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0.  ABSTRACT  (continued) 


For  the  non-linear  systems,  we  find  that  the  zero-order 
approximation  of  the  slow  variables  is  inadequate  because  it  does 
not  account  for  a non-negligible  fast  component  which  is  present 
in  these  variables.  Thus  we  develop  a new  and  simple  method  to 
correct  the  non-linear  zero-order  approximations. 
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This  thesis  applies  singular  perturbation  techniques  to  linear 
and  non-linear  models  of  a single  machine-inf inite  bus  power  system  and  a 
three  machine  power  system.  In  the  linear  realm,  we  give  a method  for 
obtaining  state  and  eigenvalue  approximations  by  computing  approximations 
to  the  block  diagonal  system  which  isolates  the  fast  and  slow  dynamics. 

We  also  give  a "growth  of  model"  method  for  determining  the  slow  and  fast 
variables  in  the  power  system  models. 

For  the  non-linear  systems,  we  find  that  the  zero-order  approxi- 
mation of  the  slow  variables  is  inadequate  because  it  does  not  account  for 
a non-negligible  fast  component  which  is  present  in  these  variables.  Thus 
we  develop  a new  and  simple  method  to  correct  the  non-linear  zero-order 
approximations . 
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1.  INTRODUCTION 

I. 1 Problem  Description 

This  thesis  is  concerned  with  the  application  of  singular  pertur- 
bation techniques  to  a class  of  power  system  dynamics.  Generally,  power 
system  models  are  of  high  dimensionality  and  are  stiff.  Therefore,  the 
motivation  for  this  study  is  clear:  since  singular  perturbation  techniques 
are  particularly  suitable  for  analyzing  stiff  systems  of  differential 
equations,  we  are  naturally  inclined  to  ask  whether  these  methodologies 
are  applicable  to  power  system  models.  This  thesis  demonstrates  that  the 
model  describing  rotor  angle  oscillations  and  associated  phenomena  can  be 
analyzed  by  singular  perturbation  techniques. 

The  singular  perturbation  method  is  often  called  a model  reduction 
method  and  it  is  in  the  sense  that  information  (e.g.  eigenvalues,  time 
responses)  about  the  system  is  obtained  from  lower  order  subsystems. 

Reduced  order  modelling  has  been  of  interest  to  power  system  engineers  for 
many  years.  Much  of  the  work  on  reduced  order  models  has  been  aimed  at 
obtaining  dynamic  equivalents  for  transient  stability  studies. 

One  of  the  "classical"  methods  of  equivalencing  [1]  merely  combines 
machine  inertias  and  parallels  their  transient  reactances.  This  method 
requires  engineering  judgment  to  determine  which  machines  can  be  represented 
by  an  equivalent  machine.  The  studies  in  which  this  method  is  often  used 
do  not  include  additional  dynamics  such  as  voltage  regulators  and  speed 
governors.  For  those  studies  in  which  the  aforementioned  dynamics  are 
important,  other  methods  have  been  devised.  One  technique  (which  can  be 
regarded  as  "classical")  is  to  neglect  "small"  time  constants  and  other 
"unimportant"  dynamics.  Another  method  [2],  which  is  limited  to  linear 
systems,  specifies  a certain  structure  for  the  reduced  model  and  certain 
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eigenvalues  of  the  full  model  which  are  to  be  retained  in  the  reduced  model. 
Then  the  gains  and  time  constants  of  the  reduced  model  are  calculated. 

This  method  requires  much  computation  and  a solution  does  not  always  exist. 

In  [3]  and  [4]  the  authors  linearize  that  portion  of  the  system  which  is  to 
be  replaced  by  an  equivalent.  Then  the  linearized  equations  are  trans- 
formed to  a diagonal  system  and  certain  "fast"  modes  and  "unimportant" 
modes  are  neglected.  This  method  requires  some  prior  knowledge  of  the 
system  in  order  to  determine  which  portion  of  the  system  can  be  replaced  by 
an  equivalent.  Also  there  is  some  judgment  involved  in  selecting  which 
modes  in  the  equivalent  are  to  be  retained  and  which  are  to  be  neglected. 
Recently,  equivalents  based  on  coherency  have  been  investigated  [5].  This 
is  a promising  avenue  of  research;  although,  again,  some  prior  knowledge  of 
the  system  is  needed  to  determine  which  machines  behave  coherently. 

Unlike  the  methods  discussed  above,  the  singular  perturbation 
approach  does  not  necessarily  generate  dynamic  equivalents.  However,  our 
approach  does  produce  reduced  order  models  in  the  sense  that  calculations 
(of,  for  example,  system  time  response)  are  performed  on  lower  order  sub- 
systems. Thus  the  full  system  is  not  evaluated  directly.  A major  advantage 
that  the  singular  perturbation  approach  has  over  other  reduced  order  modelling 
schemes  is  that  the  complete  system  response  is  recoverable  from  the  sub- 
system responses.  Thus  we  do  not  lose  information  about  the  system  as  is 
the  case  if  some  dynamics  are  completely  neglected.  Another  advantage  of 
singular  perturbation  methods  is  that  they  are  not  necessarily  restricted  to 
linear  systems. 
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1.2  Contributions 

Although  singular  perturbation  techniques  were  applied  to  a 
simple  power  system  example  in  [6],  this  thesis  appears  to  be  the  first 
application  of  these  methodologies  to  detailed  power  system  models.  A 
necessary  part  of  this  study  is  the  development  of  an  iterative  scheme  to 
obtain  approximations  to  the  block  diagonal  system  which  isolates  the  fast 
and  slow  subsystems  in  a linear  singularly  perturbed  system.  In  addition, 
we  present  a technique  (the  so  called  "growth  of  model"  method)  for 
partitioning  the  power  system  model  into  fast  and  slow  subsystems. 

Finally  we  give  a new  and  simple  method  for  correcting  the  zero-order 
approximation  of  a non-linear  singularly  perturbed  system  in  the  case  when 
the  slow  variables  contain  a non-negligible  fast  component.  All  of  these 
methodologies  are  applied  to  both  a seventh  order  single  machine- inf inite 
bus  system  and  a twentieth  order  three  machine  system. 

1.3  Review  of  Power  System  Dynamics 

Power  system  transients  range  from  travelling  wave  phenomena  on 
transmission  lines  to  very  slow  boiler  response  behavior.  In  this  thesis 
we  will  be  concerned  with  synchronous  machine  rotor  angle  oscillation 
transients  and  associated  phenomena.  These  transients  result  from  changes 
in  the  system  which  cause  a momentary  mismatch  between  generation  and  load. 
The  problem  of  sustained  mismatch  between  generation  and  load,  although  not 
treated  here,  has  received  attention  recently  [7]. 

There  are  basically  two  types  of  motion  which  can  result  from 
momentary  generation- load  mismatch,  namely  1)  stable  motion  in./(*hich  the 
generator  speeds  and  torque  angles  move  toward  the  equilibrium  corresponding 
to  the  new  system  conditions  or  2)  unstable  motion  during  which  one  or  more 
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of  the  generators  loses  synchronism  with  the  others.  The  latter  situation 
may  result  in  the  removal  of  generators  from  the  system  which  in  turn  could 
lead  to  partial  or  total  system  collapse  due  to  cascading  events. 

In  order  to  design  and  operate  the  system  in  a reasonably  reli- 
able way,  the  nature  of  power  system  dynamic  response  must  be  well  under- 
stood. This  fact  was  recognized  many  years  ago  [1]  and  the  problem  of  power 
system  "transient"  response  has  received  attention  ever  since. 
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This  stiffness  leads  to  serious  difficulties  in  carrying  out  numerical 
integration.  The  singular  perturbation  approach  alleviates  both  the 
dimensionality  and  stiffness  problems  by  analyzing  the  fast  and  slow  sub- 
systems separately.  Hence,  singular  perturbation  methods  appear  to  be 
particularly  suitable  for  studying  power  system  models. 

1.5  Chapter  Preview 

In  Chapter  2 we  discuss  some  basic  principles  of  singular 
perturbations.  We  then  present  a method  for  obtaining  an  approximation 
to  the  block  diagonal  system  which  isolates  the  fast  and  slew  dynamics  in 
a linear  system.  Finally  we  discuss  the  problem  of  partitioning  an 
arbitrary  system  into  its  fast  and  slow  subsystems  and  we  give  a procedure 
for  partitioning  our  power  system  model. 

Chapter  3 considers  a single  machine-infinite  bus  system.  First 
the  partitioning  of  this  system  is  determined  by  the  growth  of  model 
technique  described  in  Chapter  2,  Next  the  complete  linear  model  is 
analyzed  by  the  iterative  scheme  introduced  in  Chapter  2 and  both  eigenvalue 
and  time  response  approximations  are  found  to  be  quite  good.  Finally  the 
nature  of  the  zero-order  approximation  of  the  non-linear  system  is  explored. 
This  approximation  is  found  to  exhibit  some  deficiences.  A method  to  correct 
them  is  proposed  and  verified. 

In  Chapter  4 we  repeat  the  analysis  of  Chapter  3 for  a multi- 
machine system.  We  find  that  by  introducing  a change  in  angle-speed 
variables  and  by  applying  the  experience  gained  with  the  single  machine 
system,  the  partitioning  of  this  system  is  straight-forward.  State  and 
eigenvalue  approximations  of  the  linear  model  are  seen  to  be  quite  accurate 
but  more  iterations  are  required  to  achieve  this  accuracy  than  in  the  single 
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machine  case.  Finally  the  twentieth  order  non-linear  system  is  analyzed 
and  the  zero-order  approximation  is  found  to  suffer  from  the  same  type  of 
difficulty  as  the  zero-order  approximation  of  the  single  machine  system. 
The  correction  procedure  is  applied  and  the  approximations  are  seen  to  be 
improved  relative  to  the  zero-order  approximation. 
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2.  SINGULAR  PERTURBATION  METHODS 


2.1  Introduction 


Singular  perturbation  methods  [10]  are  applicable  to  systems  of 


differential  equations  of  the  following  form 
x = f (x jZ,u) , x(0)  = xQ 
M-z  = g(x,z,u),  z (0)  = zn 


(2.1a) 

(2.1b) 


where  x is  an  n-vector,  z is  an  m-vector,  u is  the  input  vector  of  dimension 
r and  H is  a small  positive  parameter  known  as  the  singular  perturbation 
parameter.  The  presence  of  M>  indicates  that  there  are  fast  dynamics  in  x 
and  z.  Settings  = 0 in  (2.1b)  is  equivalent  to  neglecting  the  fast  phenomena 
in  the  response  and  gives  the  following  slow  model 


x = f(x,z,u),  x(0)  = xc 
0 = g(x,z,u) 


(2.2a) 

(2.2b) 


where  the  overbar  indicates  the  slow  part  of  these  quantities.  The  fast  part 
of  the  response  can  be  approximately  recovered  by  solving  the  system 


Hz  = g(x,  z + z,  u + u),  z (0)  = zn  - z(0)  . 


(2.3) 


In  (2.3)  the  tilde  denotes  the  fast  part  of  these  quantities.  Notice  that 
the  approximate  responses,  x and  z + z are  obtained  by  solving  lower  order 
subsystems  thus  reducing  the  number  of  differential  equations  which  have  to 
be  integrated  simultaneously.  In  addition,  the  fast  and  slow  subsystems  are 
solved  separately  in  their  own  time  scales,  thereby  ameliorating  the  problem 
of  stiffness. 


The  linear  time -invariant  version  of  (2.1  ) is 
x = Ax  + Bz  + Gu,  x(0)  ■ xQ 
M«2  = Cx  + Dz  + Hu,  z(0)  ■ Zq 


(2.4a) 


(2.4b) 


8 

In  this  case,  the  usual  slow  subsystem  model  is 

x = (A-BD‘LC)x  + (G-BD_1H)u,  x(0)  = xQ  (2.5a) 

z = -D_1Cx  - D_1Hu  , (2.5b) 

and  the  fast  subsystem  model  is 

(iz  = Dz  + Hu,  z(0)  = Zq  - z"(0)  . (2.6) 

In  either  the  linear  or  non-linear  version  described  above,  the 
responses  x and  z are  approximated  by  x and  "z+z.  Therefore,  this  approxi- 
mation does  not  account  for  a possible  fast  component  in  x.  In  many  cases 
this  approximation  is  satisfactory.  However,  if  the  system  is  severely 
disturbed,  the  fast  response  in  x may  be  large.  In  this  case,  the  effect  of 
the  fast  variables  on  the  slow  subsystem  may  not  be  negligible.  Then  it  is 
generally  necessary  to  use  modified  procedures  to  calculate  the  approximate 
responses.  In  the  linear  case,  one  such  method  is  to  obtain  an  approximation 
to  the  block  diagonal  system  which  isolates  the  fast  and  slow  subsystems  by 
using  the  Chang  transformation  [11].  Another  method  to  arrive  at  the  same 
result  will  be  described  in  the  next  section.  In  the  non-linear  case,  the 
method  used  depends  on  the  specific  non-linear  functions  in  (2.1  ).  In 
Chapters  3 and  4 we  give  one  method  for  calculating  corrections  to  the  basic 
approximations  for  the  non-linear  equations  encountered  in  the  power  system 
models  presented  there. 

Finally,  the  singular  perturbation  parameter  p.  is  not  always  identi- 
fiable, nor  is  it  necessary  that  it  be  isolated.  The  only  requirement  for 
using  singular  perturbation  techniques  to  analyze  a system  is  that  the  system 
must  possess  the  multiple  time  scale  property.  Once  this  property  has  been 
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established  for  the  system,  then  it  is  necessary  that  the  system  be 
partitioned  correctly.  We  will  return  to  this  topic  later  in  this 
chapter. 

2.2  An  Iterative  Scheme 

One  way  to  obtain  accurate  fast  and  slow  approximations  is  to 
carry  out  the  calculations  on  a block  diagonal  system  which  isolates  the 
fast  and  slow  subsystems.  We  now  give  a method  to  construct,  approxi- 
mately, this  block  diagonal  system.  Consider  the  linear,  time- invariant 
system  (2.4  ) with  u = 0 for  simplicity 


Ax  + Bz,  x(0)  = x 


Hz  * Cx  + Dz,  z(0)  = z 


(2.7a; 

(2.7b) 


Define 


.-1, 


T)  = z + D Cx  . 


(2.8) 


We  substitute  (2.8)  into  (2.7a)  to  give 


.-1 


(A-BD  C)x  + BTL 


= A^  + BT^  . 

Now  we  take  the  time  derivative  on  both  sides  of  (2.8)  and  substitute 
(2.7b)  and  (2.9)  to  yield 


(2.9) 


* -1  -1 
M»T}1  = HD  C AjX  + (D  + HD  CB)!^ 


- C^x  + . 

This  procedure  can  be  repeated  iteratively.  At  the  end  of  the 
k-th  iteration,  the  result  is 


(2.10) 


x = Akx  + BT1k 


(2.11a) 


H^k  = 


(2.11b) 


where 


\ = Vi  • BDk-ick-i 


(2.12a) 


Ck  * “Sc-l  St- A 


Dk  = Dk.L  + P-D^  Ck_LB  . 


(2.12b) 


(2.12c) 


In  order  to  avoid  computing  ^ at  each  iteration,  we  can  approxi- 


mate these  inverses.  Let  G * Gq  + p.G^.  Then 

G_  1 = (Gq  + PG^'1 

1 2 

= Mq  + HMj.  + J ^ + * * * 


(2.13) 


It  is  clear  that  Mq  = Gq  . To  obtain  the  higher  order  terms,  we  employ 


repeated  differentiation.  For  M^  we  have 


M1  = dp,  <G0  + M'Gl^ 


(2.14) 


But  for  any  non-singular  matrix  A,  AA  =1.  Hence 


M -1  A dA_ 

dp,  dpi 


_A~^  — a”1 

dP- 


(2.15) 


Applying  (2.15)  to  (2.14)  we  obtain  for  M. 


M1  * 'go1gico1 


(2.16) 


Repeating  this  procedure  for  other  higher  order  terms,  we  find 
that  the  expansion  for  G ^ can  be  written  in  the  following  form 


g"1  = U^G^q  (I-...)]}  . 


(2.17) 


For  example,  an  approximation  of  is 


d^1  - d‘1{i-pd'lcbd'1[i-p.d‘1cbd"1(i-p.d’icbd'L)]}  . 


(2.18) 
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Notice  that  the  expression  for  contains  |i  as  a multiplying 
factor.  Hence,  at  the  end  of  each  iteration,  the  elements  of  should 
be  smaller  than  the  elements  of  Thus,  after  some  convenient  number 

of  iterations,  say  k,  the  coupling  term  C^x  can  be  dropped.  Then  the 


fast  subsystem  is  defined  as 
HT|  = fill  . 

Now  we  write  (2.11a)  as 
x = <7n  + BT| 

To  decouple  the  slow  subsystem  from  the  fast  subsystem,  we  substitute 

• 1 * 

T]  » T|  into  (2.20).  The  result  can  be  written  as 


x - M.BJ&  T|  =*  £7x 


which  suggests  the  definition  of  a new  slow  variable 


(2.19) 


(2.20) 


(2.21) 


§1  - x - M-BlS"  Tl 


With  this  definition,  (2.21)  becomes 


(2.22) 


i1  » + ^'1ti 

= ctt>l  + . 

This  scheme  can  be  employed  iteratively.  After  the  k-th 


(2.23) 


iteration,  the  result  is 


where 


5k  * ^k  + V1 


Bk  = ^k-A 


(2.24) 


(2.25) 


As  in  the  case  of  C^,  the  elements  of  are  reduced  at  each  iteration. 
Hence,  after  some  suitable  number  of  iterations,  the  coupling  terms  can 
be  dropped  and  the  slow  subsystem  defined  as 


5 - at 


(2.26) 


— ♦ 


The  number  of  iterations  performed  on  the  T1  subsystem  does  not  necessarily 
have  to  be  equal  to  the  number  performed  on  the  * subsystem.  However,  for 
computational  symmetry,  we  will  always  carry  out  the  same  number  of  iterations 
on  both  subsystems. 

We  now  turn  our  attention  to  the  problem  of  recovering  x and  z 
from  i and  "[).  First,  we  have 


\ " Vi  + Dk-ick-ix  • 


(2.27) 


If  (2.27)  is  summed  from  1 to  k,  the  result  is 


J Vl  + Vm-lVl* 

m= 1 m=l  m* 1 


(2.28) 


\-’>o  + V«-1C»-1X  ■ 

m=l 


(2.29) 


k -1 

^m-lCm-l  * 

m«l 


(2.30) 


Then  (2.29)  can  be  rewritten  as 


T1,  = z + Lx  . 


(2.31) 


In  (2.31)  we  have  dropped  the  subscript  on  T1  to  signify  that  the  result  of 
the  last  iteration  is  taken  to  be  the  true  fast  response. 

Similarly,  for  5 we  have 


’ Vl  ' “Vl*'  71 


(2.32) 


■ 1 m=l 


Z l . - M.  Z Bm  ,S'l 
, ni”  l 1 m- 1 


(2.33) 


5k-50-“  S'-I*' 

m*l 


(2.34) 
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H 


k 

I B 


m*l 

Thus,  (2.34)  becomes 


m-1 


& 


-1 


x - HHT]  , 


(2.35) 


(2.36) 


By  substituting  (2.31)  into  (2.36)  and  simplifying  the  result,  we 
arrive  at  the  following  transformation  of  variables 


I-UHL 


-HH 


(2.37) 


Incidently,  this  transformation  gives  us  a convenient  way  to  compute  the 
initial  conditions  for  simulating  the  § and  T]  subsystems.  The  inverse 
transformation  is 

I M-H 


-L 


HJ-LH 


9 

«nO 

1 

T] 

(2.38) 


This  transformation  is  a special  case  of  the  transformation  used  by 
Chang  [11];  it  is  also  obtained  in  [12]. 

The  approximation  procedure  which  has  been  outlined  in  this  section 
is  somewhat  different  than  the  basic  approximation  scheme  presented  in  the 
previous  section.  Here,  we  compute  approximations  to  the  block  diagonal 
system  which  isolates  the  fast  and  slow  subsystems.  Then,  the  original 
variables  are  reconstructed  by  means  of  the  transformation  (2.38).  The 
calculations  are  always  performed  on  lower  order  subsystems. 


i 


I 


In  later  chapters  we  use  a power  system  model  which  is  valid  for 
studying  synchronous  machine  behavior  for  a period  of  several  seconds 
following  an  initiating  disturbance.  This  model  consists  of  swing 
equations,  flux  decay  dynamics,  and  voltage  regulator  representation.  It 
possesses  the  multiple  time  scale  property  since  it  encompasses  dynamics 
ranging  from  fast  control  loops  to  system  frequency  drift  (in  the  multi- 
machine case).  Therefore,  singular  perturbation  techniques  are  appropriate 
for  studying  this  power  system  model.  By  means  of  singular  perturbations, 
we  can  obtain  reduced  order  models  which  alleviate  the  stiffness  and 
dimensionality  problems  of  the  full  order  models. 


2.4  Techniques  of  Partitioning  the  System 

One  weakness  of  the  singular  perturbation  approach  is  that  there 
is  no  general  method  for  locating  the  fast  and  slow  variables  in  an  arbitrary 
system.  Fortunately  the  model  used  in  this  thesis  allows  us  to  determine 
the  partitioning  fairly  simply.  We  shall  call  the  technique  about  to  be 
described  the  "growth  of  model"  method. 

First,  all  the  synchronous  machines  in  the  system  are  represented 
by  only  their  swing  equations  and  the  eigenvalues  of  the  linearized  system 
are  computed.  Next,  flux  decay  dynamics  are  included  with  the  swing 
equations.  The  eigenvalues  of  the  linearized  system  are  again  computed. 

It  is  an  interesting  property  of  these  first  two  models  that  the  flux  modes 
are  weakly  coupled  to  the  swing  modes.  Hence,  in  the  flux  decay  model,  it 
is  possible  to  immediately  recognize  the  flux  modes  and  the  swing  modes. 

This  behavior  is  not  coincidental  since  the  swing  modes  are  determined 
primarily  by  the  machines'  moments  of  inertia  and  the  admittances  of  the 


15 


>1 

I 


i 


interconnecting  network.  The  inclusion  of  flux  dynamics  introduces  an 
energy  dissipation  mechanism  which  somewhat  affects  the  damping  of  the 
swing  modes. 

Finally  the  voltage  regulators  are  included  with  the  model  and 
the  eigenvalues  of  the  linearized  system  are  examined.  The  swing  modes 
and  the  mode  due  to  quadrature  axis  flux  dynamics  are  discernable.  The 
other  eigenvalues  produced  by  this  model  result  from  some  interactions 
among  the  variables.  The  origins  of  these  eigenvalues  are  clarified  after 
these  interactions  are  understood. 

Thus  by  the  growth  of  model  method  we  can  identify  which  variables 
are  primarily  fast  and  which  are  primarily  slow.  This  allows  the  system  to 
be  partitioned  into  the  form  of  (2.1)  or  (2.4). 

There  is  another  method  for  partitioning  a system  when  the  system 
can  be  conveniently  represented  by  a block  diagram.  The  block  diagram  is 
examined  for  small  time  constants  and  loops  containing  large  gains  since 
their  presence  gives  rise  to  fast  dynamics.  Hence,  the  fast  and  slow 
variables  can  often  be  identified  by  inspection  of  the  block  diagram.  We 
will  use  this  method  in  the  next  chapter  as  an  aid  in  classifying  the 
dynamics  of  the  voltage  regulator  model. 
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3.  SINGLE  MACHINE  CASE  STUDY 

3 . 1 Introduction 

In  this  chapter  we  analyze  a single  machine-infinite  bus  system 
by  singular  perturbation  methods.  Both  linear  and  non-linear  results  are 
presented.  Although  this  system  has  no  true  counterpart  in  practice,  many 
practical  systems  can  be  analyzed  by  approximating  them  as  single  machine  - 
infinite  bus  systems.  In  addition  the  multi-machine  system  considered  in 
the  next  chapter  retains  many  of  the  characteristics  of  the  single  machine 
system.  Hence,  it  is  beneficial  to  analyze  and  understand  the  simpler 
system  before  trying  to  analyze  the  multi-machine  system. 


3.2  System  Description 

The  single  machine-infinite  bus  system  is  shown  in  Fig.  3.1. 

The  voltage  regulator  model  used  in  this  study  is  the  IEEE  Type  1 repre- 
sentation [13]  a block  diagram  of  which  is  shown  in  Fig.  3.2.  In  the 
model  used  here  saturation  non-linearity,  represented  on  the  block  diagram 
by  S , is  retained  but  limit  type  non-linearities  are  neglected. 

Following  are  the  differential  equations  for  each  of  the  three 
models.  The  swing  model  equations  are  well  known  [1].  The  full  model 
equations  are  derived  in  Appendix  A and  the  flux  decay  model  equations  are 
obtained  from  the  former  by  holding  the  field  voltage  (E^)  constant  and 
dropping  the  voltige  regulator  equations. 

Swing  Model 

6 = 377  <p-l)  (3.1a) 


^ = 2H  ' YE'V.sin  5] 


(3.1b) 
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Flux  Decay  Model 


e'  = t~{  - [1+  (x.-x')y]  e'  - (x,-x’ )YV,sin6  + E,,} 


'9  T 


q d 


= ^r-C(x  -x')Yvicos6  - [l+(x  -x')Y]  e^} 


6 = 377  (0-1) 


• 1 ^in 

° = 2H  C ' D(n-1)  -YVi(eqCos6  + e^sin6)] 


Full  Model 


1 *? 

h ■ r (-Rf  + ir  Efd> 

r r 


e^  = :p~{  (x^-x' )Y  V.cos6  - [l+(xq-x')y]  e^} 


6 = 377  (0-1) 


0 = ^ -D(0-1)  -YV.(e\:os6 +e^sin6)] 


Efd  ■7:‘-fKE+SE(E£d))  Efd  + V 


tKA(Rf-^E£d-Vtt',Ref)  - V 

A F 


(3.2a) 


(3.2b) 


(3.2c) 


(3. 2d) 


e’  = pH  - [l+(xd-x')Y]  e'  - (xd-x' )YVtsin6  + Efd}  (3.3a) 


(3.3b) 


(3.3c) 


(3.3d) 


(3.3e) 


(3 . 3 f ) 


(3.3g) 


V = [(l-x'Y^(e'2  +e’2)  + 2 (1-x'Y  ) (x'YV.  ) (elcos6  - e'sinfi)  + (x'YV.)2] 
t q a l a q l 


2,1/2 


(3.3h) 


SE(Efd>  = AsatexPt(Bsat)(Efd)] 


(3.3g) 
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The  parameters  used  in  this  study  are 


i 

j 


! 


D = 2.0  pu 
xd  = 1.2  pu 
x = 1.0  pu 

q 

x'  = 0.25  pu 
xg  = 0.25  pu 
T'  = 5.0  sec 

do 

T'  = 0.50  se 

qo 


«e 

*F 


0.06  sec 
0.5  sec 
1.0  sec 
25 

-0.0445 
0. 16 

0.001123 

0.3043 


3.3  Determination  of  Partitioning 

To  determine  the  fast  and  slow  variables  in  the  full  system,  we 
use  the  growth  of  model  and  block  diagram  methods  discussed  in  the  last 
chapter.  The  linearized  equations  for  each  of  the  models  described  in 
the  last  section  can  be  written  in  the  following  form 


Aw 


(3.4) 


In  each  case  the  terminal  conditions  about  which  the  linearization  is 
performed  are 


Vi  = 1.0  + jO  pu 
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0 377 

-0.2  -0.2 


the  eigenvalues  of  which  are 

-0. 1 + j8.68  • 

Since  this  model  accounts  for  only  the  mechanical  motion  of  the  rotor,  this 
mode  is  clearly  identifiable  as  the  rotor  oscillatory  mode.  Note  that  the 
"frequency"  of  this  mode  (1.38  Hz)  is  quite  typical  of  values  found  in 
practice. 

Flux  Decay  Model 

In  the  flux  decay  model 


w = [Ae'  Ae'  A6  Af2]  ' 

q d 

and  the  A matrix  is 

-0.58  0 -0.269  0 

0 -5.0  2.12  0 

00  0 377 

-0.141  0.142  -0.2  -0.28 

The  eigenvalues  of  this  matrix  are 
-0.887  + j 8.41 
-3.78 


(3.6) 


-0.301  . 


The  mechanical  mode  is  clearly  visible  (-0.887  + J8.41).  The  imaginary 
part  of  this  mode  has  not  been  affected  very  much  by  the  inclusion  of  flux 
decay  dynamics;  but  the  damping  of  this  mode  has  been  increased  substantially. 
The  energy  dissipation  provided  by  the  rotor  circuits  accounts  for  this 
additional  damping. 


..... 
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Flux  decay  dynamics  are  responsible  for  the  remaining  two  modes. 
The  root  at  -3.78  corresponds  to  Ae^  while  the  root  at  -0.301  corresponds 
to  Ae'.  This  determination  is  obvious  since  T'  = 0.5  sec.  and  T'  = 5.0 


The  quadrature  axis  rotor  coil  has  a smaller  time  constant  than 
the  direct  axis  rotor  coil  (i.e.  it  is  "faster").  Thus,  we  might  suspect 
that  more  of  the  mechanical  mode  damping  is  accounted  for  by  the  presence 
of  the  quadrature  axis  coil  than  by  the  presence  of  the  direct  axis  coil. 
If  we  modify  this  model's  A matrix  by  removing  the  dependence  of  aA  on 
Ae^  we  obtain  the  following  modified  A matrix 

f-0.58  0 -0.269  0 


-0.141 


-0.28 


The  eigenvalues  of  this  matrix  are 
-0.235  + j8.68 


-0.390  . 


Hence,  in  this  case,  a sizeable  portion  of  the  mechanical  mode  damping  is 
attributable  to  the  quadrature  axis  coil. 


Full  Model 


In  the  full  model 


w = [Ae'  AR,  Ae'  A6  Afi  AE,,  AV „] ' 
q f d fd  R 


(3.7) 


and  the  A matrix  is 


0 00  00  0.0838  2.0 

-173  417  -116  40.9  0 -66.7  -16.7 


the  eigenvalues  of  which  are 
-0.861  ± j8.39 
-3.93 

-0.362  + j0.558 
-8.53  + j8.22  . 


The  mechanical  mode  is  clearly  present  (-0.861  + j 8 . 39 ) as  is  the  mode  due 
to  (-3.93).  However,  the  origin  of  the  other  four  modes  is  somewhat 

obscure. 


To  understand  how  these  other  modes  arise,  we  make  use  of  the 
block  diagram  reduction  technique  alluded  to  in  Chapter  2.  We  consider  the 
voltage  regulator  model  of  Fig.  3.2.  The  block  diagram  of  the  linearized 
equations  correponding  to  the  loop  consisting  of  the  amplifier,  the  exciter 
and  the  Kp/T^,  feedback  path  can  be  drawn  in  the  form  of  Fig.  3.3.  The 
characteristic  equation  of  this  subsystem  is 


2 

s + 


*1*2*3 


- 1 


1 T 
L 2 


(3.8) 


r*-4 1*4 


K3 

— — — - ■ 1 - ■ 

Fig.  3.3.  Equivalent  block  diagram  of  inner  loop  of 
voltage  regulator  model. 


' 
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The  gain  is  usually  large  for  two  reasons.  First,  the  voltage 
regulator  is  a proportional  controller;  hence  large  gain  is  needed  to  produce 
small  steady  state  error.  Second,  large  gain  makes  the  regulator  faster. 

The  effect  of  increasing  gain  is  to  move  the  roots  of  (3.8)  along  the 
asymptotes  shown  in  Fig.  3.4.  Hence,  we  can  consider  this  subsystem  as 
singularly  perturbed  in  the  sense  of  high  frequency  oscillations  [6].  For 
the  numerical  values  in  the  system  matrix,  the  roots  are  -8.29  + j 7 . 95 . 

Thus  the  fast  pair  (-8.53  + j8.22)  is  attributable  to  the  aforementioned 
loop. 

To  understand  the  origin  of  the  remaining  complex  pair,  we  con- 
sider the  simpler  system  of  an  open  circuited  synchronous  machine  with 
voltage  regulator.  If  the  inner  voltage  regulator  loop  is  reduced  by 
singular  perturbation,  the  resulting  slow  subsystem  model  is  as  shown  in 
Fig.  3.5.  The  characteristic  equation  of  this  subsystem  is 

s2  + 1.45  s + 1.26  (3.9) 

the  roots  of  which  are  -0.725  + j0.857.  Thus,  we  can  conclude  that  in  the 
full  model,  the  complex  pair  -0.362  + j0.558  arises  from  the  interaction  of 
Ae'  and  ARr. 

q f 

Although  the  partitioning  of  the  full  system  into  fast  and  slow 
subsystems  is  not  unique,  one  possibility  is  to  place  -0.362  + j0.558  in 
the  slow  subsystem  and  all  other  modes  in  the  fast  subsystem.  From  the 
previous  discussions,  we  see  that  we  should  place  Ae^  and  AR^  in  the  slow 
subsystem  and  the  other  variables  in  the  fast  subsystem. 


VMI 
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Fig.  3.4.  Effect  of  increasing  gain  on  the  roots 
of  Eq.  (3.8). 


Fig.  3.5. 


Slow  subsystem  model  of  open-ci 
machine  with  voltage  regulator. 


I 
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3.4  Analysis  of  Linearized  Svstem 

In  this  section  we  analyze  the  seventh  order  linear  system  by 
means  of  the  iterative  singular  perturbation  scheme  presented  in  Chapter  2. 
First  we  examine  the  eigenvalue  approximations  for  the  seventh  order  system 
matrix  given  previously.  Table  1 shows  approximate  and  exact  eigenvalues. 
The  approximate  eigenvalues  are  computed  from  the  subsystem  matrices  after 
the  specified  number  of  iterations. 


Table  1.  Eigenvalue  Approximations  for  Seventh  Order  System 


Note  that  the  usual  zero-order  approximation  is  not  used  here 
since  it  is  not  consistent  with  the  iterative  scheme. 

After  the  first  iteration,  the  fast  eigenvalues  are  approximated 
quite  well  but  the  slow  subsystem  eigenvalues  are  somewhat  inaccurate. 
After  the  second  iteration,  all  of  the  eigenvalues  are  approximated  very 
closely. 


A very  important  test  of  the  subsystem  approximation  scheme  is 
the  accuracy  of  the  approximate  responses,  x(t)  and  z(t).  To  test  this 
aspect  of  the  methodology,  we  simulate  a system  disturbance.  The 
disturbance  chosen  for  this  study  is  a 5 cycle,  3 phase  stub  fault  applied 
at  the  infinite  bus.  Figures  3.6  - 3.12  show  the  simulation  results  for 
a period  of  10  seconds  following  the  disturbance.  In  each  case  the  solid 
curve  is  the  exact  solution  and  the  dotted  curve  is  the  approximate 
solution  obtained  by  using  one  iteration  of  the  iterative  scheme.  The 
initial  conditions  were  calculated  by  integrating  the  non-linear  equations 
(3.3)  during  the  fault  and  subtracting  from  the  variables'  values  at  fault 
clearing  the  corresponding  equilibrium  values. 

Figures  3.6  and  3.7  are  the  slow  variables.  The  predominant 
slow  behavior  of  these  responses  is  apparent.  However,  in  both  cases,  an 
initial  fast  component  is  clearly  present.  This  fast  component  is  probably 
due  to  the  forcing  effect  of  the  field  voltage  (note  that  the  field  voltage, 
E^,  appears  in  the  equations  for  e^  (3.3a)  and  (3.3b))  which  attains  a 
relatively  large  value  during  the  fault  and  shortly  after  the  fault  is  cleared. 
Once  the  field  voltage  has  decreased  to  nearly  its  equilibrium  value,  the 
slow  variables  regain  their  dominant  slow  behavior.  The  approximation  method 
used  here  takes  into  account  a possible  fast  component  in  the  slow  variables 
since  x * 5 + nHT].  This  relationship  also  shows  that  a large  initial 
condition  in  the  fast  variables  has  a non-negligible  effect  on  the  slow 
subsystem  initial  condition  because  5(0)  * x(0)  -p<HT](0). 

Even  though  there  is  some  error  in  the  approximation  of  the  slow 
variables,  it  is  clear  that  with  only  one  iteration,  the  slow  variables  are 
approximated  quite  well.  The  curves  of  the  exact  and  approximate  fast 
variables,  Figs.  3.8-3.12,  are  in  general  indistinguishable. 


0*01 


Plots 
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Since  there  are  some  discrepancies  between  the  approximate  and 
exact  responses  using  the  results  of  one  iteration,  the  system  was  simu- 
lated again  using  the  results  of  two  iterations.  However,  in  this  case  all 
of  the  approximate  responses  are  indistinguishable  from  the  exact  responses. 

3.5  Zero-Order  Analysis  of  the  Non-Linear  System 

In  this  section  we  examine  the  usual  zero-order  approximation 
(2.2)  as  applied  to  the  non-linear  system  (3.3).  The  calculations  involved 
in  the  zero-order  approximation  are  quite  simple.  We  have  previously  identi- 
fied e'  and  R,  as  slow  variables  and  e'  6,  Q,  E,., , and  V„  as  fast  variables, 

q r u rd  R 

We  assume  that  this  behavior  carries  over  to  the  non-linear  system.  Then, 
the  slow  subsystem  is 


^T“{  - [1  + (xd-x'  )y]  e'  - (xd-x'  )YVisin6  +Efd) 

do 

(3.10a) 

1 - Kp  - 

tf  ( " Rf  + tf  Efd} 

(3.10b) 

r7~{  (x  -x')YV.cos6  - [1+  (x  -x')Y]e'} 

H 1 4 d 

q0 

(3.10c) 

377  (fi  - 1) 

(3. lOd) 

1 P- 

2jj  ["H11  - YV.(e^cos6+e^sin6)  - D (H- 1 ) ] 

(3. lOe) 

tt’IKE  + SE<E£d)l 

E 

(3 . lOf ) 
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°-^!KA‘V^Efd-VV£>-'?RI 


(3. lOg) 


Vt  = [(l-x'Y)2(e’ + e^2)  + 2 (1-x'  Y)  (x'YV.)  (e^cos^  -I'sinS)  + (x'yV^2]  1/2 . 


(3. lOh) 

Equations  (3. 10c)- (3. lOh)  are  non-linear,  algebraic  relationships 
between  the  slow  variables  and  the  slow  part  of  the  fast  variables.  During 
numerical  integration  of  the  slow  subsystem  equations,  updated  values  of 
¥'  and  R^  are  used  in  (3.10c)  - (3.10h)  to  solve  for  new  values  of  e^,  6 , H, 


E,, , and  V_.  The  updated  values  of  6 and  E,,  are  substituted  into  (3.10a) 
ta  K rd 


and  (3.10b)  to  continue  the  integration  for  e'  and  Rr. 

q f 


Once  the  slow  quantities  are  known,  the  fast  subsystem  is  integrated. 
Expressing  each  fast  variable  as  the  sum  of  its  slow  and  fast  parts,  the 
differential  equations  of  the  fast  subsystem  are 


^d  = 5T-((xq-x,)YVicos(6  +6)-[l+  (x  -x')Y]  (ej  + ej)} 


(3.11a) 


5 = 377  0 


(3.11b) 


. P. 

1 r in 


H ^ — — - YV.[e'cos(5+6)  + (l-'  + e ' )sin(5  + 6) ] - DC] 

2H  Q+Q  1 * d d 


(3.11c) 


Eld-iT  t-<'t  + SEff£d+ifd»<EH+i£d>+VR  + »Rl 


(3. ILd) 


\ ■ - £ <E£d  + E£d>-’'ct’W  • vV 


(3. lie) 


V - i(l-x'Y)2  fe'2  + (el  +e')2]  + 2 (1-x’Y  ) (x'YV. ) [(el  +■  el)cos(^  + 5 ) 
t q a a i a a ’ 


- e'sin(5  +5)]+ (x’YV.)2} 1/2 


(3. Ilf) 


Since  the  slow  subsystem  has  already  been  solved,  the  slow 
variables  appearing  in  (3.11)  are  known  during  integration  of  the  fast 
subsystem.  The  fast  subsystem  is  integrated  with  a smaller  step  size  than 
the  slow  subsystem;  hence,  in  general,  it  is  necessary  to  interpolate 
between  the  known  values  of  the  slow  variables  to  obtain  the  correct  values 
to  use  in  (3.11). 

To  explore  the  nature  of  this  zero-order  singular  perturbation 
analysis  of  the  non-linear  system,  we  simulate  the  system  response  for  the 
same  disturbance  used  in  the  last  section.  Figures  3.13-3.19  show  the 
results  of  this  simulation.  In  each  case,  the  solid  curve  is  the  exact 
response  and  the  dashed  curve  is  the  zero-order  approximate  response.  It  is 
quite  clear  that  these  responses  are  very  similar  to  their  linear  counter- 


parts. 


We  have  previously  noted  that  the  slow  variables  contain  a fast 


component.  Since  the  zero-order  approximation  of  the  slow  variables  does 

not  admit  a fast  component,  e'  and  R,  are  not  approximated  very  well  by  ¥' 

q f q 

and  R^.  Despite  this  inaccuracy,  the  fast  variables  are  approximated  quite 
well.  Thus,  we  are  motivated  to  search  for  a procedure  to  correct  the  slow 
variable  approximation. 

3 .6  A Correction  Method  for  the  Non-Linear  System 

One  method  to  obtain  corrections  for  non-linear  systems  is  given  in 
[14].  Another,  simpler  method  to  approximately  account  for  a fast  component 


in  the  slow  variables  will  now  be  given.  Consider  the  system 


x = f(x,z),  x(0)  = xQ  . 


(3.12) 


I 


ONM MM 


Fig.  3.19.  Plots  of  exact  and  zero-order  approximation  of  V 
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Suppose  that  z is  written  as 
z = cp(x)  + z. 


(3.13) 


where,  symbolically,  cp(x)  is  the  solution  of  g(x,z)  = 0 and  zf  is  the 


fast  part  of  z.  Substituting  (3.13)  into  (3.12)  we  get 


x = f(x,cp(x)  + zf)  . 


(3.14) 


Let  f be  such  that  it  can  be  decomposed  as  follows 

f(x,cp(x)  + zf ) = f1(x,cp(x))  +f2(cp(x),  zf ) 

= f L (x)  + f2(cp(x),  z f ) . (3.15) 

As  an  approximation,  let  us  use  a previously  obtained  solution  for  x and 


z,  in  f . Let  these  solutions  be  denoted  by  x°  and  z . Then  (3.15) 


becomes 


x = fL(x)  + f2(cp(x°),zf°) 


(3.16) 


The  integral  equation  form  of 

(3.16)  is 

t 

t 

x-xQ  = J fL(x)dr  + 

J f2(cp(x°),  Zf°)dT. 

(3.17) 

0 

0 

Let 

k = J f,  (cp(X°),  2 °)dT. 


(3.18) 


We  write  (3.17)  as 


= XQ  + k+  J f1(x)dT  + J f £ (cp (x°)  , zf°)dT  - k 
0 0 


(3.19) 


1 

A 
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As  an  approximation,  neglect  the  fast  part  of  x in  the  first  integral.  Then 
the  slow  part  of  x is  given  by 

xs  = fL(xs),  xg(0)  = xQ  + k , (3.20) 

and  the  fast  part  of  x is  approximated  by 

t 

xf  - -k  + J f2(cp(x°),  zf°)  dT  . (3.21) 

0 

Notice  in  (3.20)  how  the  presence  of  fast  variables  in  the  slow  equations  can 
change  the  slow  initial  conditions. 

Once  the  slow  subsystem  has  been  solved  with  the  corrected  initial 
conditions,  the  fast  subsystem  can  be  resolved  using  the  corrected  values  of 
x,  namely 


Lzf  = g(xg+xf,  zf  + zg),  zf  (0)  = z0"zs  (0) 


(3.22a) 


0 = g(xs,zs)  . 


(3.22b) 


Note  that  this  scheme  could  be  employed  iteratively. 

Applying  this  method  first  to  (3.3a)  we  find  that  the  slow  part  of 
e'  is  obtained  as  the  solution  of 

q 

e'  = pK  - [1+  (xd-x')Y]  e'  - (x^x'  )YV.sin6s  + Efd  }, 


e'  (0)  = e' (0)  + k. 

qs  q i 

and  the  fast  part  of  e'  is 

q 

(x  -x')YV.  t 

e'  = -k — i j [sin6  (-l+cos6f  ) + 

q,  i i , s z 


*0  0 
o . . o. 


cos6  sin6c  ]dT  + —7-  j E , , dr 

s t 1,0  rd, 

d0  0 f 


(3.23) 


(3.24) 


49 


and 


(x.-x')YV,  » 


kl  = • 


T,  — ur  [sin6s°(-l  + cos6f°)  + 


cos6s°  sinfif°]dT  + j E°d  dT. 

d0  0 f 


(3.25) 


Similarly  for  (3.3b),  the  slow  part  of  R^  is  solved  from 


Rf  ■ £ (-RE  + F|E£d  >•  R£  (0)  ■ Rf<°>  + k2 
s F s F s s 


(3.26) 


and  the  fast  part  is 


Rff  = _k2  + ~2  I Efdf  dT 
f TF  0 f 


(3.27) 


where 


k2  ' “2  I Efd_  dT  ' 


(3.28) 


V 0 ‘ f 


6g  and  E^d  appearing  in  (3.23)  and  (3.26)  are  solved  from  the 
s 


same  algebraic  equations  as  in  the  original  approximation.  The  integrals 
in  (3.24),  (3.25),  (3.27),  and  (3.28)  can  be  evaluated  by  any  convenient 
numerical  integration  scheme  since  the  arguments  appearing  under  the  integral 
sign  are  known  functions  of  time.  We  also  note  that  it  is  not  necessary 
to  carry  out  the  integrals  in  (3.25)  and  (3.28)  to  infinite  upper  limits. 


From  Fig.  3.16  and  3.18  we  observe  that  6^  and  E°d  both  become  approximately 


zero  for  t > 5 seconds. 

Finally  the  fast  subsystem  is  solved  as  in  (3.11)  but  this  time, 


both  the  slow  and  fast  parts  of  e^  and  are  substituted  instead  of  just 


the  slow  part. 

To  test  this  procedure,  we  use  the  same  system  disturbance  as  in 
previous  sections.  In  carrying  out  the  indicated  calculations,  the  zero- 
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order  approximation  is  obtained  first.  The  zero-order  fast  solution  is 
used  to  evaluate  the  necessary  integrals.  Using  the  corrected  initial 
conditions,  the  slow  subsystem  is  solved  again.  To  complete  the  pro- 
cedure, the  corrected  responses  e'  + e'  and  R.  + R,  are  substituted 

q q,  f f, 

Ms  n s f 

into  the  fast  subsystem  equations  and  this  subsystem  is  solved  again. 

The  results  of  the  simulation  are  shown  in  Figs.  3.20-3.26.  In 
general  the  agreement  between  the  exact  and  approximate  responses  is  good. 
Regarding  the  slow  variables,  it  is  clear  that  the  correction  procedure 
has  not  only  accounted  for  the  fast  component  but  also  corrected  the 
initial  conditions.  Using  these  corrected  slow  solutions  during  inte- 
gration of  the  fast  subsystem  yields  generally  excellent  fast  approximations. 
Hence,  in  this  case,  the  correction  procedure  for  non-linear  systems 
performs  adequately. 

3 .7  Conclusions 

In  this  chapter  we  have  verified  that  the  single  machine-infinite 
bus  system  possesses  the  separation  of  time  scales  property.  Using  the 
growth  of  model  and  block  diagram  techniques,  we  have  discovered  which 
dynamics  of  the  model  are  primarily  responsible  for  each  mode.  It  was 
found  that  e^,  6,  H,  Ec^,  and  constitute  a suitable  fast  subsystem 
while  e^  and  R^  comprise  the  slow  subsystem. 

We  have  applied  linear  singular  perturbation  methods  to  the 


linearized  seventh  order  model.  In  particular,  one  iteration  of  the 
iterative  technique  produces  quite  acceptable  eigenvalue  and  time  response 
approximations.  Two  iterations  give  almost  perfect  approximations.  In  the 


Fig.  3.21.  Plots  of  exact 


Fig.  3.24.  Plots  of  exact  and  corrected  approx' 
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simulations,  we  noted  that  the  slow  variables  can  have  a significant  fast 
component.  The  iterative  scheme  for  linear  systems  automatically  accounts 
for  fast  components  in  the  slow  variables. 

Finally,  we  investigated  the  behavior  of  the  zero-order  singular 
perturbation  approximation  for  the  non-linear  model.  We  found  that  the 
forcing  effect  of  the  field  voltage  can  cause  discrepancies  between  the 
true  and  approximate  slow  responses.  A method  to  approximately  account  for 
the  inclusion  of  fast  dynamics  in  the  slow  subsystem  was  proposed.  This 
technique  proved  successful  for  the  particular  system  disturbance  considered 
here. 
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4.  STUDY  OF  MULTI-MACHINE  SYSTEM 

4.1.  Introduction 

In  this  chapter  we  generalize  the  results  of  Chapter  3 by  studying 
a multi-machine  system.  We  find  this  system  to  be  very  similar  to  the 
single  machine  system.  In  particular,  the  variables  which  we  previously 
identified  as  slow  and  fast  remain  in  these  respective  subsystems.  But  we 
also  find  an  additional  slow  mode  due  to  the  system  frequency  drift.  The 
properties  of  the  angle-speed  modes  are  clarified  when  the  original  angles 
and  speeds  are  transformed  to  some  new  angles  and  speeds. 

We  present  both  linear  and  non-linear  results.  The  linear  itera- 
tion technique  does  not  encounter  any  difficulties  with  this  system.  How- 
ever, for  the  particular  disturbance  considered,  the  slow  variables  generally 
contain  a significant  fast  component.  Thus,  the  zero-order  non-linear 
approximation  suffers  from  the  same  difficulty  as  in  the  single  machine 
case.  The  method  presented  in  Chapter  3 to  correct  the  slow  approximation 
is  found  to  improve  the  approximations  relative  to  the  zero-order  results. 

4.2.  System  Description 

The  system  studied  in  this  chapter  is  the  9 bus,  3 machine  system 
shown  in  Fig.  4.1.  The  transmission  line  data  are  given  in  Table  2 and  the 
synchronous  machine  parameters  are  given  in  Table  3.  The  exciter  parameters 
are  the  same  for  all  three  machines  and  are  identical  to  the  settings  used 
in  Chapter  3.  The  disturbance  used  throughout  this  chapter  is  a 3 phase,  j 

5 cycle  fault  applied  near  bus  8 which  is  cleared  by  opening  the  line  from 
bus  8 to  bus  9. 


j 


2 3 

1.63  ou 

(D 

0.0654 ou—  ~ 
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Fig.  4.1.  3 machine,  9 bus  cesc  syscem. 

Note:  All  load  flow  information  is  for  predisturbance  conditions. 
The  base  power  is  100  MVA. 
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Table  2.  Transmission  System  Data 
(100  MVA  Base) 


Line  # 

From 

To 

R(pu) 

X(pu) 

B/2  (pu) 

1 

1 

4 

0 

0.0567 

0 

2 

4 

5 

0.017 

0.092 

0.079 

3 

5 

6 

0.039 

0.170 

0.179 

4 

3 

6 

0 

0.0586 

0 

5 

6 

7 

0.0119 

0.1008 

0.1045 

6 

7 

8 

0.0085 

0.072 

0.0745 

7 

8 

2 

0 

0.0625 

0 

8 

8 

9 

0.032 

0.161 

0.153 

9 

9 

4 

0.01 

0.085 

0.088 
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Table  3.  Synchronous  Machine  Data 
(100  MVA  Base) 


Parameter 

Machine  # 

1 

2 

3 

xd (pu) 

0.6 

0.8958 

0.9 

xq (pu) 

0.58 

0.8645 

0.85 

x^(pu) 

0.056 

0.110 

0.18 

x^ (pu) 

0.0608 

0.1198 

0.1813 

X 

►o 

C 

— ** 

0.0608 

0.1198 

0.1813 

q 

T'  (sec) 
d0 

4.0 

6.0 

5.0 

T’  (sec) 

0.25 

0.54 

0.65 

^0 

ra(pu) 

0 

0 

0 

H(sec) 

23.64 

6.4 

3.01 

D(pu) 

9.6 

2.5 

1.0 
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We  now  discuss  the  choice  of  angle  and  speed  variables  for 
describing  this  system.  The  well  known  non-linear  equations  for  the  swing 
model  are 


Q. 

i 


377(0.-0^, 


2H . ^Pin. 

i l 


D. 


(^.-1) 


3 

- Z 

j-1 


i-2,3 


V ! V'.Y , .cos  (0 
i j ij 


ij 


+ 6.1-  6tl)J  , 


i - 1,2,3  . 


(4.1a) 


(4.1b) 


When  (4.1  ) is  linearized  about  the  equilibrium  point  corresponding  to  the 
post-disturbance  network,  the  resulting  system  matrix  is  found  to  be 


0.0 

0.0 

-377.0 

377.0 

0.0 

0.0 

0.0 

-377.0 

0.0 

377.0 

0.0144 

0.0235 

-0.203 

0.0 

0.0 

-0.131 

0.0930 

0.0 

-0.195 

0.0 

0.216 

-0.371 

0.0 

0.0 

-0.161 

The  eigenvalues  of  this  matrix  are 
-0.0857  + j 12 . 90 
-0.0978  + j6 . 09 
-0.197  . 

The  mode  -0.197  is  the  well  known  frequency  drift  mode  which  corresponds  to 
the  common  movement  of  all  of  the  machine  speeds  to  the  new  equilibrium 
speed.  The  two  complex  pairs  are  rotor  oscillation  modes.  Since  one 
"frequency"  is  more  than  twice  the  other  "frequency",  we  suspect  that  these 


|j 

J 


modes  may  be  attributed  to  something  other  than  individual  rotor  oscillations. 
In  particular,  with  the  removal  of  the  line  from  bus  8 to  bus  9,  machine  2 
is  relatively  weakly  connected  to  machine  1.  However,  the  coupling  between 
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machines  2 and  3 remains  relatively  strong-  Thus,  it  is  likely  that  one 
of  the  oscillatory  modes  represents  the  common  oscillation  of  machines  2 
and  3 with  respect  to  machine  1;  while  the  other  mode  is  the  oscillation 
of  machine  2 with  respect  to  machine  3.  This  situation  is  similar  to  the 
one  described  in  [6]. 

To  more  clearly  exhibit  the  behavior  described  above,  we  intro- 
duce the  following  change  of  variables.  First,  since  the  system  frequency 
drift  involves  all  of  the  rotating  inertia  of  the  system,  we  define  the 
system  frequency  to  be 


where 


0 

r 


H1Q1 


*1*2 


H3n3 


(4.2a) 


H = Hl  + H2  + H3  . 


(4.2b) 


Second,  to  display  the  combined  motion  of  machines  2 and  3,  we  define 

6 and  Cl  as 
c c 


5c  = hH  (H2621  + H3531) 


Cl  = 

C 


H2n2  + hn3 


Q 


where 


h = H2  + H3  . 


(4.3a) 

(4.3b) 

(4.3b) 


Finally,  the  oscillation 
by 


of  machine  2 with  respect 


to  machine  3 is  described 

(4.4a) 

(4.4b) 


These  variables  are  similar  to  the  ones  used  by  Stanton  in  [15]. 
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In  these  variables,  the  system  matrix  is 


An 


A6 


AQ 


A6 


AO 


d •- 


-0.198 

0.00776 

0.00486 

0.000733 

-0.00181 

0.0 

0.0 

377.0 

1 

1 

0.0 

0.0 

0.0122 

-0.113 

-0.191 

» 

1 

-L 

0.0304 

-0.00454 

0.0 

0.0 

0.0 

• 

« 

0.0 

377.0 

_ -0.0292 

0.164 

-0.0292 

i 

• 

-0.426 

-0.175 

The  eigenvalues  of  the  indicated  diagonal  blocks  are,  respectively  -0.198, 
-0.0954  + j6.54,  and  -0.0877  + jl2.6S.  These  values  correspond  quite 
favorably  with  the  eigenvalues  of  the  original  matrix  indicating  that  this 
choice  of  variables  is  more  appropriate  for  describing  the  system  than  the 
original  set  of  angle-speed  variables. 

Using  the  new  angle-speed  variables  along  with  the  basic  equations 
derived  in  Appendix  A,  the  complete  set  of  equations  for  this  system  is 


P. 

in. 


P. 

in. 


P. 

in. 


Qr  " S'' 


n - n 

r c 


n +n  + -r1  n, 

r c h d 


n +n  - n. 

r c h d 


R H 

0(Pr-i)  + (Dx  i - d2-d3)Oc  + (-D2  -if  + 03  i )"d 


■ E (e'  i + e’  i ) 1 
i-1  qi  qi  di  di 


V 


Rf  " T (_Rf.+  T Efd.}  » 
1 Fi  1 Fi  1 


i - 1,2,3 


qi 


Td 


Oi 


(xd.’  Xl}id  + Efd.] 

ni  i i i 


i - 1,2,3 


(4.5a) 


(4.5b) 


(4.5c) 
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where 


D = D1  + d2  + d3 


(4.5k) 


2 2 

V = v,  + v 
C.  , d.  q. 

1 T i i 


(4.51) 


v,  = e'  + x'  i 
di  di  1 qi 


(4.5m) 


v = e ' - x!  i, 

a . q . l a . 
4i  4i  i 


The  currents  appearing  in  (4.5  ) are 


(4.5n) 


V Yn<'i1cos9Li  • eilsln6u) 


+ 'f12led2OOS(?12+57Sc+^  V • 

- eLsinlf> u+  £ V it 


(4.6a) 


i<,1  = Yu(ei1cos9n  + eilsl”9u) 

+ Yi2l8i2cos<eu+  6C+  iT  V + ed2Sin<?12+  V tT  8d)] 

+ Yute'  cos<eu+ A«c-  6d)  + T 6d)l 


(4.6b) 


% ‘ Y22(8i2tOS  S22  - e;2Sln622) 


4 V21(a'iCOs<e21-  ^ Sc-  ^ 5d)  - •;i*inP2l-  ^ V TT  V1 


4 V32[a^cos(932-  sd)  - 6d>] 


(4.6c) 
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i = Y (e’  cos  9 + e'  sin  0 ) 

q2  22  q2  22  d?  22 

H H 

+ Y21[eq1COS(32r  5c'  T V + ed1SLn(e21*  5c'  IT  5d)] 


+ Y23[e>  cos(023-  5d)  + sin(323-  5d)] 


(4.6d) 


% = Y33(ed3COS  933‘ 

H H 

+ hl[ed]cos®3l-  57  6c+  if  V - e;isin<331-  h7  V if  '5d)] 


+ Y32[e^cos(032+  6d)  - e^sin(032+  6d)] 


(4.6e) 


i = Y..(e'  cos  9,,  + e'  sin  3,,) 
q9  33  q3  3_>  d3  33 


+ Y31[e^cos(931-  f-  5c+  6d)  + e^sin(0 


31-S76c+lf 


Y32[ei2COS(8  32+  5d}  + ed9Sin(332+  5d)]  ’ 


(4 . 6 f ) 


4.3  Determination  of  Partitioning 

In  the  previous  section,  using  the  swing  model,  we  found  the 
location  of  the  swing  and  drift  modes.  The  next  step  in  the  growth  of 
model  method  for  determining  the  system  partitioning  is  to  examine  the 
eigenvalues  of  the  flux  decay  model.  Unfortunately,  the  flux  decay  model 
does  not  have  an  equilibrium  for  the  post-disturbance  network  chosen  for 
this  study.  Thus,  the  linearized  equations  cannot  be  obtained.  However, 
the  complete  model  (4.5)  does  possess  an  equilibrium  and  the  eigenvalues  of 
the  linearized  equations  fall  into  groups  which  are  reminiscent  of  the 
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location  of  eigenvalues  in  the  single  machine  case.  Hence,  we  can  use  our 

experience  with  the  single  machine  system  to  determine  which  dynamics  are 

primarily  responsible  for  each  mode.  Having  made  this  determination,  we 

can  proceed  to  partition  the  system  into  slow  and  fast  subsystems. 

The  system  matrix  for  the  complete  model  is  shown  in  Fig.  4.2. 

The  eigenvalues  of  this  matrix  are 

-1.26  + j 12 . 95 
-8.16  + j 7 . 6 9 
-8.56  + j8.26 
-8.50  + j8.09 
-0.988  + j 5 . 35 
-6.90 
-4.45 
-1.90 

-1.22  + jO. 989 
-0.245 

-0.505  ± jO. 76 1 
-0.289  + J0.480  . 

Based  on  our  examination  of  the  swing  model  in  the  last  section,  it  is 

clear  that  the  pair  -1.26  + j 12 . 95  is  primarily  due  to  46^  and  A while 

the  pair  -0.988  + j5.35  arises  from  A6,  and  Aft,.  The  real  root  -0.245  is 
— ad 

the  system  frequency  drift  mode  and,  hence,  is  mainly  determined  by  Afi  . 

Now,  using  the  experience  gained  during  our  study  of  the  single 
machine  system,  it  is  clear  that  the  pairs  -8.16  + j7.69,  -8.56  + j 8 . 26 , 
and  -8.50  + j8.09  are  mainly  produced  by  the  interaction  of  AE,,  and  AV 

IQ  R 

in  the  three  individual  voltage  regulators.  Furthermore,  the  pairs 
-1.22  + j0.989,  -0.505  + j0.761,  and  -0.289  + j0.480  obviously  arise  from 
the  combination  of  Ae^  and  AR^  in  the  respective  machine  models.  Finally, 
the  three  Ae^j 1 s are  primarily  responsible  for  the  three  real  roots,  -6.90, 
-4.45,  and  -1.90. 

▼ 

I 


System  matrix  for  full  model. 

Note:  only  non-zero  entries  shown. 
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It  is  apparent  that  this  system  behaves,  in  many  ways,  like 

three  single  machine  systems.  Thus,  for  partitioning  the  system  into 

fast  and  slow  subsystems,  the  variables  should  be  grouped  per  phenomenon 

(as  opposed  to  per  machine)  using  the  previously  identified  fast  and  slow 

variables  as  a basis  for  the  grouping.  Of  course  we  must  include  the  new 

variable,  A£)  , which  does  not  appear  in  the  single  machine  system. 

Based  on  this  discussion  and  the  system  eigenvalues,  we  choose 

as  the  slow  subsystem  the  variables  Afi  , AR  , AR  , AR  , Ae ' , Ae'  , and 

r fl  f2  f3  ql  q2 


The  remaining  variables  constitute  the  fast  subsystem. 


4.4  Linear  System  Study 

The  eigenvalue  approximations  of  the  full  system  along  with  their 
exact  counterparts  are  shown  in  Table  4.  We  notice  that  most  eigenvalues 
are  approximated  very  closely  after  the  second  iteration.  However,  the 
approximations  corresponding  to  exact  modes  at  -1.22  + j0.989  and  -1.90  are 
not  completely  satisfactory  after  only  two  iterations.  In  fact,  it  takes 
six  iterations  to  obtain  reasonable  approximations  for  these  modes. 

Regarding  this  latter  behavior,  we  observe  that  -1.22  + j0.989 
is  the  largest  of  the  slow  eigenvalues  while  -1.90  is  the  smallest  of  the 
fast  eigenvalues.  The  ratio  of  their  magnitudes  is  1.21;  hence,  there  is 
no  clear  separation  between  them.  Based  on  this  fact,  we  can  speculate 
that  these  eigenvalues  are  somewhat  sensitive  and  therefore  more  difficult 
to  approximate  than  the  others.  In  particular,  these  eigenvalues  converge 
to  their  correct  values  non-uniformly . 

The  next  eigenvalue  to  the  left  is  -4.45.  The  separation  between 
-4.45  and  -1.90  is  more  favorable  than  the  separation  between  -1.90  and 
-1.22  + j0.989.  Hence  we  are  led  to  believe  that  if  the  variable  which  is 


l 


Table  A.  Exact  and  Approximate  Eigenvalues  for  20th  Order  System 


responsible  for  the  -1.90  mode  were  placed  in  the  slow  subsystem,  this 


non-uniformity  of  convergence  might  not  occur.  By  examining  the  Ae^  block 


in  the  system  matrix,  we  see  that  the  eigenvalue  -1.90  is  primarily  deter- 


mined by  Ae'  . So  a more  reasonable  choice  for  the  slow  subsystem  would 
d3 


include  Ae^  with  the  variables  already  selected. 


The  simulations  reveal  some  interesting  characteristics  of  this 


system.  The  time  responses  for  the  post-disturbance  network  are  given  in 


Figs.  4.3  - 4.22.  The  initial  conditions  for  these  linear  simulations  were 


calculated  by  the  same  procedure  as  was  used  to  calculate  the  linear  initial 


conditions  in  Chapter  3.  In  each  case,  the  solid  curve  is  the  exact  solution 


while  the  dotted  curve  is  the  approximation  obtained  using  the  results  of 


three  iterations  of  the  iterative  procedure.  It  was  found  that  if  four  or 


more  iterations  were  performed,  the  exact  and  approximate  curves  were 


virtually  indistinguishable. 


We  note  that,  as  in  the  single  machine  case,  the  slow  variables 


(Figs.  4. 3-4. 9)  contain  fast  components.  In  addition,  it  appears  that  the 


machine  1 slow  variables  contain  more  of  the  rotor  angle  oscillations 


(particularly  the  lower  frequency  oscillations)  than  their  machine  2 and  3 


counterparts.  This  observation  suggests  that  the  mode  -1.22  + j0.989  is 


a machine  1 mode.  Thus,  Ae'  and  ARf  respond  more  strongly  to  the  fast 

ql  fl 


phenomena  than  the  machine  2 and  3 slow  variables.  The  effect  of  large 


field  voltage  (^E^)  is  readily  observed  in  the  Ae^'s  and  AR^'s  as  was  the 


case  in  the  single  machine  system.  The  slow  part  of  all  of  the  slow 


variables  is  clearly  present. 


The  response  of  A fir  (Fig.  4.3)  is  very  interesting.  A Cl  is 


proportional  to  the  mismatch  between  total  system  real  power  generation 


and  total  load.  Approximately  the  first  4 seconds  of  the  A H response  is 


Fig.  4.5.  Plots  of  exact  and  approximate  AR 


Plots  of  exact  and  approximate  Ae’ 


Fig.  4.11.  Plots  of  exact  and  approximate  Ae' 
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dominated  by  fast  dynamics.  This  behavior  can  be  explained  as  follows. 

In  this  study,  the  loads  are  modelled  by  constant  impedances.  Thus,  the 
real  power  absorbed  by  the  system  is  very  sensitive  to  the  internal  machine 
voltages  and  angles.  Since  both  the  internal  voltages  and  angles  have  fast 
parts,  the  total  system  load  should  have  a fast  part;  hence,  we  should  not 
be  surprised  that  the  system  frequency  has  a fast  part.  In  addition,  once 
the  internal  voltages  and  angles  have  come  nearly  to  steady  state  (after 
about  4 seconds),  then  Aftf  is  dominated  by  a single  exponential  characteristic 
(namely  the  pole  at  -0.245).  Once  the  voltages  and  angles  reach  steady 
state  conditions,  then  the  power  absorbed  by  the  system  also  comes  to 
steady  state.  Thereafter,  Aft^  is  governed  primarily  by  the  system  speed 
regulation  characteristic. 

The  fast  subsystem  responses  also  have  some  interesting  features. 

Ae'  (Fig.  4.10)  contains  a much  larger  component  of  the  rotor  angle 
dl 

oscillations  than  do  Ae'  and  Ae'  . Thus,  we  are  led  to  believe  that  the 

d2  d3 

mode  -6.90  is  primarily  due  to  Ae^  . Incidently,  the  other  machine  1 fast 
variables  also  contain  more  rotor  angle  oscillation  component  than  their 
counterparts  of  machines  2 and  3. 

Regarding  the  rotor  angle  oscillations,  A6c  and  Aflc  (Figs.  4.19 
and  4.20)  appear  to  be  very  pure,  namely  they  consist  primarily  of  the  lower 
frequency  oscillations.  A6d  and  Aft^  (Figs.  4.21  and  4.22),  on  the  other 
hand,  show  evidence  of  a mixture  of  the  two  oscillatory  modes.  However, 
the  lower  frequency  mode  appears  to  dominate  indicating  that  the  inter- 
machine oscillation  is  not  excited  very  much  by  this  disturbance. 
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4.5  Study  of  the  Non-Linear  System 

In  this  section  we  consider  the  non-linear  system  (4.5).  Based 
on  the  linear  study  of  the  previous  section,  we  choose  ftr,  , and  e^ 

(i * 1,2,3)  as  the  slow  variables  and  the  remaining  states  as  the  fast 
variables.  The  calculation  of  the  zero-order  approximation  proceeds  along 
the  same  lines  as  in  the  single  machine  case.  To  test  the  zero-order 
approximation,  the  system  was  simulated  for  a period  of  10  seconds 
following  the  test  disturbance.  The  results  of  the  simulation  are  shown 
in  Figs.  4.23  - 4.42.  As  before,  the  solid  curve  is  the  exact  solution  and 
the  dashed  curve  is  the  approximation. 

We  notice  the  similarity  between  these  responses  and  the  corre- 
sponding responses  of  the  linearized  system.  Thus,  the  general  comments 
made  about  the  physical  interpretation  of  the  linear  responses  can  be  carried 


over  to  the  non-linear  system. 

The  zero-order  approximation  of  the  slow  variables  (Figs.  4.23  - 
4.29)  evidently  suffers  from  the  same  difficulties  as  the  zero-order 
approximation  of  the  single  machine  slow  variables,  namely,  absence  of  a 
shift  of  the  slow  initial  conditions  and  a lack  of  fast  components.  The 
approximations  of  the  fast  variables  also  have  some  shortcomings.  In 
particular,  in  several  cases  it  is  apparent  that  the  slow  part  of  the  variable 
is  in  error;  furthermore  the  "frequency"  of  the  slower  of  two  rotor  angle 
oscillation  modes  is  not  correct.  Since  this  mode  is  a component  of  many 
other  responses,  this  error  propagates  through  most  of  the  zero-order  fast 
approximations . 

Thus  we  again  find  it  necessary  to  compute  corrections  to  the  zero- 
order  approximation.  To  this  end,  we  will  use  the  method  of  Chapter  3. 

Most  of  the  labor  involved  in  applying  this  method  to  the  slow  subsystem 


zero-order  approximation 


zero-order  approximation 


exact  and  zero-order  approximation 


zero-order  approximat 


zero-order  approximation  of 


Fig.  4.29.  Plots  of  exact  and  zero-order  approximation  of 


zero-order  approximation  of 


exact  and  zero-order  approximation  of 


Fig.  4.33.  Plots  of  exact  and  zero-order  approximation  of  V 


Plots  of  exact  and  zero-order  approximation  of  E 


zero-order  approximation  of 


zero-order  approximati 


Plots  of  exact 
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equations  is  contributed  by  the  machine  currents.  The  expressions  for 
the  currents  are  given  below  in  terms  of  the  original  angular  variables. 

This  is  done  for  notational  and  computational  simplicity.  It  is  an  easy 
matter  to  calculate  the  slow  and  fast  parts  of  the  original  angles  in  terms 
of  the  respective  parts  of  6^  and  6^  since  these  two  groups  of  angles  are 
linearly  related. 

The  components  of  the  currents  are 
3 

i , = E Y.  . [e ' cos  (9  .+6.1-6  )-el  sin©  . . + 6 -5..)]  (4.7a) 

d±  j = 1 iJ  d.  ij  jl  xV  ij  jl  il'J 


i = Z Y . . [e ' cos  (9..+6.1-5.1)+e'  sin  (6. .+6. .-6.,)] 
j=1  iJ  qj  J-J  Jl  il  d.  ij  jl  il'J 

We  write  e'  and  6.,  as  the  sum  of  their  slow  and  fast  parts 

dj  Jl 


ed.'  *1  +ed.„  - 

J js  jf 

jl  j Is  jlf’ 


j = 1,2,3 


j =2,3 


(4.7b 


(4.8a) 


(4.8b) 


We  substitute  (4.8)  into  (4.7)  and,  after  some  manipulation,  obtain  the 
slow  and  fast  parts  of  the  currents 


i,  = L Y . . [e ' cos  (8  +6..  -6.,  ) - e'  sin(9..+5  *6.,  )]  (4.9a) 

d is  j=1  IJ  djs  lj  Jls  Lls  qj  Ij  Jls  ils" 


J 

i,  = 2 Y . , { [e ' cos  (9.. +6..  -6..  )-e'  sin©  +6..  -6.,  )]  X 

dif  j=l  ij  djs  1J  jls  lls  qj  LJ  Jls  lls 


[-l  + cos(fl  -6  )]  - [(e-  + ed  )sin(9ij  +6jls  ‘ 6ils) 

J js  j£ 


+ e^cos(9..  +6  jIs  - 6ils)]sin(6jlf  - 6.lf) 

+ e’ _fCos(9..  +6  jls  - 6tls)cos(6jlf  - 6ilf)} 


(4.9b) 


I 


■ 
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1 = E Y.  . [e  cos  (9  . . + 6 . . -6.,  ) + e sin(9  . . + 6 . , -6.,  )] 

q.  . , ij  1 q.  v 11  jls  lls'  d.  ij  jls  ils'J 

M1S  J = 1 J j j s J J 


(4.9c) 


i = E Y..{[e'  cos(9  +6..  *6.,  ) 
qif  j-1  LJ  qj  ij  jls  Lls 


+ e^jsSLn^i-i  + 6Jls  " 6ils)]  t'L  + cos(6jif  " 6ilf)] 

+ [(e’  + ed.  >COs(9ij+6jls-6ils)-ei.Sin(9ij+6jls*6ils)]sin(6jlf-6ilf) 

J s J i J 


+ e'  sin(9..+6.1  -6..  )cos  (6  . . , - 6 . . ,) } 

' ij  jls  ils  jlf  llf  J 


Approximation  of  Equation 

In  order  to  apply  the  method  of  Chapter  3 to  (4.3a),  we  first 
approximate  the  quotients  as 


(4 . 9d ) 


in. 


0 -hfi  /H. 
r cl 


P,  P. 

in.  . in. 

— L_  + Jl  — I n 

n H.  n 2 c 

r l n 

r 


(4.10a) 


' in„ 


(1  +STi  /h 

r c 3d 


Pin.  Pin. 

T <V113tVh> 


n 


' in- 


r Q 


' in. 


(4.10b) 


’ inn 


n +0  - H-O./h 
r c 2d 


IT”  TT^c-Wh> 


n 


(4.10c) 


1 


Equation  (4.10)  follows  since,  generally,  and  are  small.  Now  we 

write  e'  , i,  , and  i as  the  sum  of  their  slow  and  fast  parts.  The  sums 
di  di  qi 

in  (4.5a)  can  be  split  as  follows 


■■■■■■■ 


i 
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slow:  E (e ' i + e'  i , ) 

1-1  qi  qis  dis  dis 


(4.11a) 


fast:  E [e1  i + (e'  + e'  )i,  + e'  i,  ] . 

i-1  qi  qif  dis  dif  dif  dif  dis 


(4.11b) 


Finally,  since  0^^  = Q = 0,  the  slow  part  of  0^  is  approximated  by  the 
solution  of 

T>  J.  P J.  P 

3 


P.  + P,  + P. 

• 1 ^*^7  in« 

°™  - 25  ( 07 1 - Dpr,-l>  <•:  /,  + >1- 


rs 


i=l  Mis  Mis  is  is 


Qrs^0)  = + C1  +c2  +c3  +c4  + c5  + c6  +C7 


and  the  fast  part  is  approximated  as 


rf 


P.  „ n° 
h Lnl  Qcf 

*ci  + ir^rJ  -^2"dT  - c 
1 0 “r 


2 2H 


Pin.  t fi°,  + H,f)°  /h 

2 I < Cf.o2  3 > « 


in3  f , °cf  ' H2ndf/h 


Q 


"3  p , cf  ~'2  df'  \ , 1 h p Oo  . 

'C3  2H  J ( Oo2  )dT  * c4  + H,  2H  J ^cfdT 


Q 


‘C5  ' 2H  J (Qcf + H3Qdf/h)dT  ' C6  ' 2H  I (ncf  ' H2Qdf/h)dT 

0 0 


i t 3 

1 f r „ r .0  .0  ✓ | O i O N , O ,0  .0  i i « 

+ V >ld.„  + V,  V 1]iT 

0 i=l  if  is  if  if  if  is 


where 


. pio.  * 0° 


•1  2H  0 „02 

u r 


Pin,  » n°,+H,n°  /h 

I (~'-;0F'  ) dT 


2H 


Q 


(4.12) 


(4.13) 


(4.14a) 


(4.14b) 


(4.14c) 


in.  * ^ - HJ1  /h 

-1 ; < cf.o22df  > * 


2H 


n 


1 h [>“0  .. 

C4  ' 5T  25  I ncf  dT 

1 0 


(4 . 14d ) 


C5  * ‘ 2H  i*  (ncf +H3ndf/h)  dT 
0 


(4 . 14e ) 


'6--™l"<Pcf-“2nd£/h>dT 
0 


(4. 14f  ) 


® 3 


7=  - t Z [e'°  i°  + (e'°  + e ' ° )i°  + e ' 0 L°  ]} 

7 2HJ0  i-1  qi  qif  dis  dif  dif  dif  dis 


dT 


(4. 14g) 


Approximation  of  Equations 


In  (4.5c)  we  write  i and  Ec.  as  the  sum  of  their  slow  and  fast 

di  fdi 


parts.  The  slow  part  of  e'  is  given  approximately  by  the  solution  of 


1 1 


’i.  Ti 


6q.  ‘ <V  Xi>1d.  + Efd.  >• 


Oi 


Hs  i 


is  is 


(4.15) 


V°>  ” + C8  "C9 


and  the  fast  part  is  approximated  by 


qif 


-C8  - 


(x  -x>) 
i 


Td 


Oi 


I C£d* 
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C9  + TI  L EfdifdT 


(4.16) 
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(xd.'  Xi> 
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1°  dr 
dif 


(4.17a) 


(4.17b) 


d0i  0 if 

Approximation  of  Equations 

Writing  E, , as  the  sum  of  its  slow  and  fast  parts  in  (4.5b) 
i 

gives  the  following  for  the  approximation  of  the  slow  part  of  Rf 

^ j 

V. 


Rf  T Rf  + t Efd  ^ ’ 

is  Fl  is  F.  is 


Rf  (0)  = Rf  (0)  + c . 

is  i 


(4.18) 


The  fast  part  is  approximated  by 

*F.  t 


Rf./  'C10  + *2  J Efd.fdT 

if  T 0 if 

Fi 


(4.19) 


where 


i (•  _o 

C10  = T„ 2 J Efd . 


dT 


(4.20) 


‘Fi"  0 if 


The  results  of  the  calculations  are  shown  in  Figs.  4.43  - 4.62. 
These  plots  were  made  after  3 passes  through  the  iterative  cycle.  Regarding 
the  slow  variables  (Figs.  4.43-4.49),  it  is  clear  that  a fast  component  has 
been  added,  but  the  "amount"  of  fast  component  is  not  enough. 

The  fast  variables  (Figs.  4.50  - 4.62)  also  show  some  improvement 


over  their  zero-order  counterparts.  In  particular,  the  "frequency"  of  the 
slower  rotor  angle  oscillatory  mode  is  more  accurate  in  the  corrected  case. 
The  increased  accuracy  of  this  mode  favorably  affects  many  of  the  other  fast 


corrected  approximation  of 


Fig.  4.46.  Plots  of  exact  and  corrected  approximation  of 


Fig.  A. 49.  Plots  of  exact  and  corrected  approximation  of 


exact  and  corrected  approximation  of 


variables.  In  addition,  it  appears  that  the  slow  parts  of  some  of  the 
fast  variables  are  more  accurate  due  to  the  improvement  of  the  slow  sub- 
system approximation. 

Thus,  it  is  obvious  that  while  this  procedure  has  produced  some 
improvement  over  the  zero-order  approximation,  the  improvement  is  not  as 
dramatic  as  in  the  single  machine  case.  The  probable  reason  for  this 
behavior  is  that  the  separation  between  the  fast  and  slow  dynamics  is  not 
as  great  in  the  multi-machine  system  as  it  is  in  the  single  machine  system. 

We  observed  this  fact  when  studying  the  linearized  system.  Hence,  approxi- 
mating the  slow  subsystem  as  an  integrator  of  the  fast  dynamics  is  probably 
not  very  accurate.  This  state-of-af fairs  might  be  improved  somewhat  if 

e'  were  included  with  the  slow  subsystem. 
d3 

4.6  Conclus ions 

In  this  chapter,  we  first  explored  eigenvalue  location  for  two 
models  of  the  three  machine  system,  namely,  1)  swing  model  and  2)  full 
twentieth  order  model.  Using  the  swing  model,  we  were  able  to  identify  the 
system  frequency  drift  mode  and  two  rotor  angle  oscillatory  modes.  We 
found  that  the  smaller  of  these  modes  was  primarily  determined  by  the 
combined  inertia  of  two  machines  while  the  larger  of  the  two  was  associated 
with  the  oscillation  of  one  of  these  machines  relative  to  the  other.  This 
observation  prompted  us  to  transform  the  angular  variables  to  more  readily 
display  this  behavior. 

The  eigenvalues  of  the  full  model  were  found  to  fall  into  groups 
the  locations  of  which  were  very  close  to  the  locations  of  some  of  the  eigen- 
values of  the  single  machine  system.  Hence  we  were  able  to  identify  which 
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phenomena  were  responsible  for  each  of  the  eigenvalue  groups.  We  con- 
cluded that  to  achieve  time  scale  separation,  it  is  necessary  to  group 
the  equations  on  a per  phenomenon  basis. 

The  variables  Q , e'  ,e'  , e'  . R_  , R,  , and  R.  were  chosen 

r ql  q2  V V f2  f 3 

to  comprise  the  slow  subsystem  with  all  other  variables  in  the  fast  sub- 
system. With  this  split,  there  is  hardly  any  separation  between  the  fastest 
of  the  slow  modes  and  the  slowest  of  the  fast  modes.  Because  of  this  lack 
of  separation,  the  eigenvalue  approximations  took  six  iterations  to  con- 
verge to  reasonable  accuracy.  The  state  approximations  showed  good 
accuracy  after  three  iterations  and  excellent  accuracy  after  four  iterations. 

As  in  the  single  machine  case,  the  slow  variables  exhibited  a 
large  component  of  fast  phenomena.  Thus  the  zero-order  non-linear  approxi- 
mation of  the  slow  subsystem  was  not  particularly  accurate;  however,  the 
fast  subsystem  approximation  was  reasonably  accurate.  Applying  three 
iterations  of  the  non-linear  correction  procedure  resulted  in  some  improve- 
ment in  the  approximations.  But  this  improvement  was  seen  to  be  less  dramatic 
than  in  the  single  machine  probably  because  of  the  relatively  small  separation 
between  the  slow  and  fast  subsystems. 


I 


143 


5.  CONCLUSION 


This  thesis  applies  singular  perturbation  techniques  to  a power 
system  model  describing  rotor  angle  swings  and  associated  phenomena. 

First  we  considered  a single  machine- infinite  bus  system.  Using 
the  linearized  model  we  ascertained  the  proper  partitioning  of  the  system 
into  fast  and  slow  subsystems.  We  also  found  that  the  eigenvalues  and 
states  of  the  linearized  model  could  be  approximated  quite  closely  by 
using  a block  diagonalization  procedure.  Finally  the  zero  order  approxi- 
mation of  the  non-linear  model  was  investigated.  The  slow  subsystem 
approximation  was  found  to  be  in  error  because  it  failed  to  account  for  a 
large  fast  component  in  the  slow  variables  due  to  forcing  action  of  the  field 
voltage.  A method  to  correct  for  this  inadequacy  was  introduced  and  was 
seen  to  perform  well. 

We  concluded  by  studying  a three  machine  system.  From  the 
linearized  swing  model,  we  found  that  this  system  possessed  a slow  mode 
which  was  not  present  in  the  single  machine  case,  namely,  the  system 
frequency  drift.  In  addition,  the  nature  of  the  swing  eigenvalues  suggested 
a change  of  angle-speed  variables  in  order  to  facilitate  partitioning  the 
system.  Once  this  change  of  variables  was  made  and  after  the  equations  were 
grouped  on  a per  phenomenon  basis,  partitioning  the  system  into  fast  and  slow 
subsystems  was  relatively  easy.  Again  we  found  the  state  and  eigenvalue 
approximations  of  the  linearized  model  to  be  very  good,  although  more 
iterations  were  required  in  this  case  because  of  the  relatively  small 
separation  between  slow  and  fast  subsystems.  Finally  the  zero-order  approxi- 
mation of  the  twentieth  order  non-linear  system  was  tested.  The  true  slow 
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subsystem  responses  were  seen  to  have  a large  fast  component  as  in  the 
single  machine  case.  We  applied  the  correction  procedure  but  with  less 
dramatic  success  than  in  the  single  machine  case  probably  because  of  the 
poor  separation  exhibited  by  this  system. 

Directions  for  future  research  are  numerous.  An  automated 
methodology  for  partitioning  an  arbitrary  system  should  be  sought.  More 
work  needs  to  be  done  regarding  correction  methods  for  non-linear  singularly 
perturbed  systems.  In  this  same  light,  models  of  some  power  system  components 
such  as  voltage  regulators  and  governors  contain  such  "hard"  non-linearities 
as  limiters  and  deadbands.  Methodologies  need  to  be  developed  to  account 
for  these  non-linearities. 
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APPENDIX  A:  Derivation  of  Synchronous  Machine  Model 

Figure  Al  shows  a simplified  three  phase  synchronous  machine 
which  illustrates  some  of  the  conventions  and  nomenclature  which  will  be 
used  in  the  following  discussions. 

Two  circuits  are  shown  on  the  rotor.  The  coil  labelled  f is  the 
main  field  winding;  the  coil  labelled  g is  a ficticious  coil  often  intro- 
duced [16]  to  account  for  so-called  "transient"  effects  in  the  quadrature 
axis.  Its  presence  is  an  attempt  to  model  eddy  currents  which  can  flow 
in  the  solid  rotors  of  high  speed  turbo-alternators  under  non-steady  state 
conditions.  Additional  rotor  circuits  are  often  included  [9]  to  account 
for  "subtransient"  effects.  These  rotor  circuits  are  needed  in  order  to 
model  eddy  currents  flowing  near  the  surface  of  the  rotor  (or  in  the 
amortisseur  winding  of  a salient  pole  machine)  which  contribute  to  the 
damping  of  rotor  oscillations.  Hence,  the  additional  rotor  circuits  are 
often  included  when  studying  those  problems  in  which  an  accurate  assessment 
of  rotor  oscillation  damping  is  considered  essential.  Since  this  thesis  is 
concerned  with  time  scale  separation  properties  and  not  with  accurate 
calculation  of  synchronous  machine  damping,  these  additional  rotor  circuits 
have  not  been  included. 

Notice  also  that  the  quadrature  axis  is  shown  in  Fig.  Al  to  be 
leading  the  direct  axis.  This  state-of-af fairs  was  the  formerly  adopted 
convention  [16];  but  some  recent  literature  [17]  has  suggested  the  opposite 
convention.  Neither  convention  appears  to  have  any  particular  advantage 
over  the  other. 

Based  on  the  discussion  of  the  previous  paragraphs,  the  following 
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per  unit  equations  for  the  machine's  electrical  variables  can  be  written  [18] 


Vd  "3£*d 


v = — 'll  + 0 ^ 
q u)q  q d 

0 = r i + — if 

g g u>0  g 


v,  = r,i,  + — ill, 
f f f f 


(Ala) 


(Alb) 


(Ale) 


(Aid) 


*d  = "Vd  + Mdf  Lf 


(Ale) 


'll  = -L  i + M i 

q q q qg  & 


(Alf) 


*f  ■ + Lt  4 


wig) 


ilf  * -M  i + L i 

g qg  q g g 


(Alh) 


The  symbolism  used  here  is  explained  in  Appendix  B. 

Some  additional  comments  concerning  these  equations  are  in  order. 
First,  the  stator  resistance  has  been  neglected  for  simplicity.  Second, 
the  zero  sequence  equations  are  not  included.  The  zero  sequence  relation- 
ships are  only  needed  when  unbalanced  operation  is  being  considered.  All 
the  studies  carried  out  in  this  thesis  are  for  balanced  operation.  Third, 
saturation  of  the  synchronous  machine  is  neglected.  It  has  been  found  [18] 
that  saturation  can  be  approximately  accounted  for  empirically.  The  usual 
procedure  is  to  adjust  the  mutual  inductances  as  a non-linear  function  of 
the  currents.  One  of  these  approximate  methods  could  have  been  employed 
in  this  thesis;  however,  because  our  main  emphasis  is  studying  time  scale 
separation,  we  chose  the  simplest  synchronous  machine  model  which  neverthe- 
less retains  the  basic  time  scale  information. 

Now  we  carry  out  some  manipulations  and  simplifications  to  put 
the  general  equations  (A1  ) in  their  final  form.  First  the  and 
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terms  in  (Ala)  and  (Alb)  are  neglected  and  Q is  set  to  Qq  in  the  same 
equations.  The  former  step  is  almost  always  taken  [19].  The  inclusion 
of  these  terms  accounts  for  certain  rapidly  decaying  components  of  the 
stator  current  during  transients.  The  torque  produced  by  these  components 
of  the  current  is  generally  negligible  except  in  the  case  of  a fault  on 
the  generator  terminals  in  which  case  the  neglected  components  of  the 
braking  torque  can  be  accounted  for  empirically.  A more  compelling  reason 
to  neglect  these  terms  is  that  if  these  terms  were  included  in  the  stator 
equations,  then  similar  terms  arising  from  inductances  in  the  network  to 
which  the  machine  is  connected  should  be  included.  For  large  networks, 
this  requirement  is  clearly  prohibitive. 

The  latter  simplification  carried  out  above  is  often  done  [19]. 

As  pointed  out  in  [20]  both  the  terminal  voltage  and  network  inductive 
reactance  vary  in  direct  proportion  to  the  machines'  frequencies.  Hence, 
in  a network  which  is  mainly  inductive,  relatively  small  variations  in 
machine  speed  do  not  appreciably  effect  the  calculation  of  network  currents. 

The  modified  stator  voltage  equations  are 


vd  = (A2a) 

vq  = noyd  . (A2b) 


Now  solve  (Alh)  for  1 and  (Alg)  for  i yielding 

§ ^ 

1 M 

i - T~  i|r  + 7^  i (A4a 

8 Lg  8 Lg  q 

1 Mdf 

4 " rf  *f + —f  Ld  • (A4b 

Finally  substitute  (A4a)  into  (A3a)  and  (A4b)  into  (A3b)  to  give 

M 

vd  ' Vi1,  ' °0  *g  <A5a 

v,  ' -"oLdld  ‘“iTJ'i  <S5b 

where  _ 

Ld  " Ld  - -tj-  <A6> 

is  the  direct  axis  transient  inductance  and 

M2 

L'  - L - -p-  (A7) 

q q L 

g 

is  the  quadrature  axis  transient  inductance. 

At  this  point  we  let 

xd  ■ Vd  . <A8a 

x’  - Qor  (A8b 

M 

ed  = 'no  +g  <A8c 

9 

•i  ■ % if  *f  • <A8d 

The  first  two  quantities  above  are  respectively  the  direct  axis  and 
quadrature  axis  reactances.  Thus  using  (A8) , (A5)  can  be  written  in  final 
form  as 
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v,  * x'  i + e ' 
d q q d 

v = -x' i , + e 1 . 

q d d q 


(A9a) 


(A9b) 


We  will  discuss  the  role  of  (A9  ) after  we  derive  the  remaining  equations 

for  this  synchronous  machine  model. 

First  we  obtain  the  differential  equations  for  e'  and  e'. 

d q 

Substituting  (A4a)  into  (Ale)  gives 

d\|r  r r 

+ -°-  (Al0) 

Multiply  this  equation  by  -Q.M  /L  . Then  by  using  the  definitions  of 

0 qg  g 

e'  and  L',  the  result  can  be  written  as 

d q 

1 ded  r r 

77  + 7s  e!.  - fl-(L  -L')i  =0  . (All 

(jo«  at  L a L 0 a a a v 

0 g g H M 4 


(All) 


The  ratio  r /L  is  the  per  unit  quadrature  axis  open  circuit  time  constant, 

O O 

T'  .We  can  also  make  the  substitution 
q0(pu) 


x -x'  = ft- (L  -L' ) 
q q 0 q q 


(A12) 


With  these  substitutions,  (All)  becomes 


_L_  1 

dt  " t; 


0 (pu) 


[ -e'  + (x  -x')i  ] . 

d q q q 


(A13) 


Finally,  we  multiply  both  sides  of  (A13)  by  uu^.  The  ratio  uon/T' 

0 0 q 


1/T'  if  T'  is  expressed  in  seconds.  Thus,  the  final  form  of  the 


0(pu) 


differential  equation  for  e'  is 

d 


ded  i 

dT  = ~[‘e'd  + VW 


(A14) 
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L 


Carrying  out  similar  manipulations  on  (Aid)  gives  the  following 


differential  equation  for 


de^ 


M. 


dt 


^T”[  - e ' - (xd-x^)idJ  + u )Qno  vf  . 

d„  f 


(A15) 


This  equation  (or  one  similar  to  it)  is  often  written  in  a slightly 
different  form  which  we  now  obtain. 

In  the  per  unit  system  being  used  in  this  thesis  [21],  base  field 
current  is  chosen  so  that  one  per  unit  field  current  establishes  the  same 
fundamental  component  of  the  mutual  flux  as  does  one  per  unit  stator 
current.  In  other  words,  the  base  field  current  is  chosen  such  that 


M 


df 


Va(b)Vf(b) 


M 


df 


(A  16) 


where  the  tilde  is  temporarily  used  to  denote  per  unit  quantities  and  the 
subscript  "b"  indicates  base  quantities.  Md^  is  the  difference  between  the 
per  unit  direct  axis  inductance  and  the  per  unit  stator  leakage  inductance. 


It  is  common  practice  [13]  to  choose  the  base  exciter  voltage  so  that 


one  per  unit  exciter  voltage  produces  one  per  unit  open  circuit  terminal 
voltage  on  the  air  gap  line  of  the  synchronous  machine.  From  (A3b)  the 
open  circuit  voltage  is 

Vq  = Wf  * (A17) 

Note  that  (A17)  is  in  ordinary  units  (volts,  amps,  etc).  Hence,  we  set 

E, 


'fd(b) 

v . * uu„M,  , L 

a(b)  0 df  r. 


(A18) 
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in  other  words 


Va(b)  = ;JJ0Mdf 
Efd(b)  rf 


(A19) 


But  from  (A16)  we  have 


0JnM  f = .,a ^ ^ M, f 

0 df  if(b)  df 

and  from  the  definition  of  per  unit  quantities 


(A20) 


rf  = rf  i 


f (b) 


(A21) 


f(b) 

Substituting  (A20)  and  (A21)  into  (A19)  gives  upon  rearrangement 


vf  (b)  = Mdf 


(A22) 


fd(b)  vf 

Now  let  the  voltage  applied  to  the  field  circuit  be  Efd (volts).  This 
value  can  be  written  as  the  product  of  its  per  unit  value  and  the  base 
exciter  voltage 

E = E E . (A23) 

^fd  rd  fd(b) 

To  put  this  quantity  in  per  unit  on  the  base  field  voltage,  we  divide 
(A23)  by  v b)  to  give 


Jfd 


f (b) 


= E 


‘ f d (b ) 


fd  v 


(A24) 


f(b) 


Substituting  (A22)  into  (A24)  and  using  T'  = Lf/r  yield 

Q0 


'fd 


= E, 


Jf  1 


Vf(b)  ^fd~Mdf~T' 


(A25) 
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In  other  words,  a per  unit  value  of  (which  might  come  from  the  solution 

of  the  exciter  equations)  must  be  multiplied  by  in  order  to  refer 

it  to  the  field  circuit. 

We  now  return  to  (A15)  and  hereafter  drop  the  tilde  convention 
for  per  unit  quantities.  We  substitute  the  per  unit  value  of  from 

(A25)  into  (A15)  to  give 


dt  = tT_[  ‘ 6q  ’ (VXd):Ld  +Efd] 
do 


(A26) 


where  all  quantities  are  in  per  unit  except  T'  which  is  in  seconds.  This 

do 

final  form  follows  since  = 1 per  unit. 

We  next  consider  the  mechanical  equation  of  motion  which  can  be 
written  as  [18] 


2H  = T.  - T , - D<p-1)  . 

dt  in  ele  ' 


(A27) 


The  last  term  on  the  right  hand  side  of  (A27)  is  a damping  term  which  is 
often  introduced  [9]  to  account  for  such  effects  as  changes  in  turbine 
power  with  speed  and  changes  in  connected  load  with  system  frequency.  The 
developed  electrical  torque  is  [18] 


T . = • ♦ i. 

ele  d q q d 


(A28) 


Now  if  (A4b)  is  combined  with  (Ale)  ijt^  can  be  written  as 


*d = -Ldid  + ~fh 


(A29) 


Similarly,  ’|r  is 


f = -L'i  + -r3^  \|( 

q q q 


(A30) 
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Having  derived  the  differential  equations  for  the  synchronous 
machine  model,  the  importance  of  the  currents,  i^  and  i , is  clearly  seen. 
Equation  (A9)  is  one  relationship  between  the  machine  terminal  voltage  and 
current.  There  is  another  constraint  between  the  voltage  and  current, 
namely  the  constraint  imposed  by  the  network  to  which  the  machine  is 
connected.  Therefore,  the  procedure  for  computing  the  currents  in  terms 

of  the  state  variables  depends  on  how  the  network  is  modelled. 

• • 

A condition  which  results  from  neglecting  the  and  if  terms 
in  (Ala)  and  (Alb)  is  that  the  network  currents  and  voltages  are  in  a quasi- 
sinusoidal  steady  state.  The  reason  we  say  "quasi-steady  state"  is  that 
the  per  unit  frequency  of  the  voltages  produced  by  the  synchronous  machines 
is  equal  to  their  per  unit  speed.  Since  these  speeds  vary  during  a transient, 
the  currents  and  voltages  cannot  truly  be  considered  as  sinusoidal.  This 
non-sinusoidal  character  is  always  neglected  in  large  system  studies. 

If  the  currents  and  voltage  are  considered  to  be  sinusoidal,  then 
they  can  be  represented  as  phasors.  The  frequency  of  these  sinusoids  is 
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taken  to  be  the  fundamental  frequency  (60  Hz).  In  general,  the  only  speed 
effect  which  might  be  included  is  the  variation  of  the  magnitude  of  the 
machine's  terminal  voltage  with  speed. 

Regarding  the  phasors,  the  formerly  accepted  practice  [16]  was 
to  take  direct  axis  quantities  as  real  quantities  and  quadrature  axis 
quantities  as  imaginary  quantities.  More  recently  [17],  the  opposite 
convention  has  been  advocated.  In  this  thesis,  we  will  use  the  former 
convention.  For  example,  the  phasor  for  the  machine's  terminal  voltage 
is  written  v^  + jv  . Note  that  phasors  written  on  this  basis  are  local  to 
a particular  machine  because  they  depend  on  the  position  of  the  machine's 
direct  axis  for  reference.  The  implication  of  this  observation  is  that 
the  angular  differences  between  the  machines ' direct  axes  are  fundamental 
in  calculating  the  network  currents.  The  variation  of  these  angular 
differences  during  transients  is  therefore  of  great  importance  in  assessing 
system  stability  following  disturbances. 

A vast  simplification  in  computing  the  machine’s  currents  is 
obtained  if  transient  saliency  is  neglected;  that  is,  we  let  x^  = x^  = x'. 
Hereafter,  we  assume  that  this  has  been  done.  Note  that  (A9)  can  now  be 
written  as  one  equation 


e^  + je^  = jx'dd  + jiq)  + vd  + jvq  . (A34) 

Equation  (A34)  is  sometimes  referred  to  as  the  voltage  behind  transient 
reactance  model.  Using  (A34)  allows  us  to  easily  compute  the  components 
of  the  currents  in  terms  of  the  state  variables  because  the  transient 
reactance  can  be  combined  with  any  external  impedance. 


! 
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As  mentioned  previously,  the  task  of  calculating  the  currents 
depends  very  much  on  how  the  network  is  modelled.  In  this  thesis,  we 
represent  the  loads  by  constant  impedance.  With  this  simplification,  the 
network  can  be  reduced  so  that  only  the  internal  generator  nodes  remain. 
Then  the  generator  currents  can  be  written  as 


I = Y V 

(1)  (1) 


(A35) 


where  1^  and  V ^ are  vectors  the  components  of  which  are  the  complex 
generator  currents  and  voltages  respectively.  The  subscript  1 indicates 
that  the  reference  for  these  phasors  is  the  direct  axis  of  machine  1.  It 
is  easy  to  show  that  any  machine  phasor,  say  a,  which  uses  the  machine  1 
direct  axis  as  reference  can  be  referred  to  its  own  direct  axis  by  applying 
the  transformation 

-js. 


jl  - 

a , = e J a 

(m) 


(A36) 


where  is  the  angle  between  the  direct  axes  of  machine  j and  machine  1 

and  the  subscript  m denotes  machine  reference. 


Now,  write 


I = T I 

(1)  L(m) 

V.. , - TV.  , 
(!•)  (m) 


(A3  7 a) 
(A3 7b) 


where 


T = 


1 + jO 
0 


j6 


21 


0 

0 


nl 


(A3  8) 
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Substituting  (A37)  into  (A35)  gives 


I,  = T 1YTV,  . 
(m)  (m) 


X*YT  V 


(m) 


Y V 


(m) 


(A39) 


where  the  superscript  * indicates  complex  conjugate.  If  we  write 


Y.  . = Y. . e 
30  30 


j0 . . 

30 


(A40) 


then  the  elements  of  Y are 


Y. . = Y. . e 
30  3.J 


j(9ij  + 6jl'6il) 


(A41) 


Splitting  (A39)  into  its  real  and  imaginary  parts  gives 


i,  + ji  = (Y  + jY.)(e’  + je')  . 
d q r i d J q 


(A42) 


In  (A42),  i^,  i^,  e^,  and  e^  are  vectors;  the  subscript  r denotes  "real 
part  of"  and  the  subscript  i denotes  "imaginary  part  of."  Equating  real 
and  imaginary  parts  on  both  sides  of  (A42)  yields 


A - A | 

i,  = Y e'  - Y.e’ 
d r d l q 


(A43a) 


i = $ e'  + Y.e'  . 
q r q id 


(A43b) 


Che  i-th  component  of  (A43a)  can  be  written  as 


2 (?  e'  - e'  ) 


d . r d 

3-  J = 1 i-J  J 


l.  . q . 

30  J 


(A44) 


which,  if  we  use  (A41),  becomes 


i , = E Y [e'  cos  (0  . + 5 . -6  . . ) - e ' sin(9  + 6 . , -6  . , ) ] . (A45) 

di  j=l  1J  dj  ij  11  qj  1J  jl  Ll 
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where  n is  the  number  of  generator  nodes.  Similarly  for  (A43b)  we  get 


i = E Y . . [e ' cos  (8  . . + 6 . , -6  . , ) + e ' sin (9  . . + 6 . , -6  . , ) 3 . 
qi  j=l  LJ  qj  1J  Jl  Ll  dj  1J  jl  11 


(A46) 


From  (A45)  and  (A46)  the  importance  of  the  angular  differences  is 
clearly  seen.  These  angles  must  be  included  in  the  state  description  of 
any  system  of  synchronous  machines.  Differentiating  the  angle  6^  with 
respect  to  time  gives 


>ji = wv 


= 377  (0.-0^,  j = 2 , 3 , • • • ,n  . 


(A47) 


The  factor  377  is  needed  to  convert  from  per  unit  to  radians.  With  (A47) 
we  complete  the  derivation  of  the  general  equations  for  the  synchronous 
machine  model . 

The  single  machine -inf inite  bus  system  shown  in  Fig.  3.1  allows 
some  simplification  of  (A45)  and  (A46 ) . Choosing  as  the  reference 


and  writing 


j (x'  +x£) 


(A48) 


it  is  easy  to  show  that 


i , = Y(e ' + V.  sin6) 
d q I 

i = Y (-e ' + V . cosS ) 
q d l 


(A49a) 


(A49b) 


The  other  equations  needed  for  both  the  single  machine  and  multi- 
machine models  are  the  equations  describing  the  voltage  regulator.  By 


inspection  of  Fig.  3.2,  these  equations  are 
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APPENDIX  B:  Notation 


With  reference  to  the  synchronous  machine  electrical  equations  of 
Appendix  A,  the  following  notation  has  been  adopted: 

Subscripts 

d - direct  axis  quantities 
q - quadrature  axis  quantities 
f - main  field  quantities 
g - quadrature  axis  "field"  quantities 
Other  Symbols 


i - current 


L - self  inductance 


M - mutual  inductance 


r - resistance 


v - voltage 

>]/  - flux  linkage 

u)q  - machine  rated  electrical  speed  in  radians/second 
Cl  - machine  per  unit  electrical  speed 
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APPENDIX  C:  Listing  of  Computer  Programs 


This  appendix  gives  the  listings  of  the  three  major  programs  used 
during  this  research. 

1.  Program  LINSP  does  the  calculations  for  a linear  singularly 
perturbed  system. 

2.  Program  VRSMSP  does  the  calculations  for  the  non-linear  model 
of  the  single  machine-infinite  bus  system. 

3.  Program  MMSP  does  the  calculations  for  the  non-linear  model 
of  the  multi-machine  system. 

These  programs  call  the  following  system  library  subroutines: 

1.  RKGS  - from  SSP  for  integrating  differential  equations. 

2.  LEQTLF  - from  IMSL  for  solving  linear  equations. 

3.  EIGRF  - from  IMSL  for  calculating  eigenvalues. 

4.  INITT 

5.  BIN ITT 

6.  NPTS 

7.  XFRM 


8.  YFRM 


CHECK 


10.  DSPLAY 


From  Tektronix  Software  Package  for  drawing  graphs. 


11.  FRAME 


12 .  ANMODE 


13.  CPLOT 


14.  LINE 


15.  DLIMY 
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2H  C 7 , 1 3 ) , wOOT  (20 ) , 0$  ( 1 3 , 1 3)  , wK  ARE  A ( 1 3 ) , AS  (7 , 7 ) , BN  (7  , 
313) ,TR (20.2a) ,y»i (2^,30) ,XI (7) ,eta( 13) ,PLT (501) ,TFST 
a(50l) ,TSL0(5tfM ,XIS (7, 5011 , ETAS fl 3,501 ).WH (501), WS (20, 
5501) , Tw (501 ) ,01 (13. 13) 

COMPLEX  UA(7) , UO ( 1 3 ) ,Z(t,  l) 

C0MM0N/8LK 1 /PA , N 

C0MM0N/8LK2/ A , NS 

C0MM0N/8UK3/0.NF 

CUMMON/Bl«a/WS, TW, NSW 

CUMM0N/BLK5/XIS, TSLO.NSS 

COMMON /5LKfc/ET AS , TFST.NSF 

EXTERNAL  FW.OW, FSLn.OSLO.FPST.OFST 

1 F0RMaT(aD 

4 FORMAT (Afl) 

u FQRMAT(I) 

31  FORMAT(G) 

40  Format ( * ??',ai,'??  Type  again  ',s) 

100  TYPE  200 

200  FORMAT ( ' GIVE  A MATRIX  FILE  NAME  '.$) 

ACCEPT  4.NNAM 
TYPE  3^0 

300  F0RMAT('  NS  a ',S) 

ACCEPT  It, NS 
TYPE  400 

400  FqrmAt('  NFs  «,S) 

ACCEPT  11, NF 
NsNS+NF 

CALL  IFILE(20,NNAM) 

00  520  I * 1 » N 

READ (20 , 3 1 ) (FA(I, J) , J = 1,N) 

500  CONTINUE 

00  1300  1 = 1 , N 

1 I * I “NS 

00  1200  J = l , N 
JJ=J-NS 

IF  C 1 1 J 603 , 608 , 900 
600  IF  (J J)  700,700,800 

700  A(I,J)=FA(I,J) 

GO  TO  1200 

800  8(I,JJ)*Fa(I,J) 

GO  TO  1200 

900  IF(JJ)  1000,1000,1100 

1000  CU (II , J) »FA  (I, J) 

GO  TO  1200 

1103  0 (I  I , JJ) =FA  ( I , J) 

1200  CONTINUE 

1300  CONTINUE 

type  1430 

1400  FORMAT!'  GIVE  IC  FILE  NAME  ',$) 


ACCEPT  4 1 NN am 
rtLL  TPTI  F f R0  . NN AM) 
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nu  15J0  I*l,N 
»£AUC2C',3ll  ^ICCI) 

1500  CONTINUE 

TYPE 

1600  FORMAT  (*  number  of  itepationSs  *,$} 

ACCEPT  It. NITER 
TYPe  173R 

1 7aa  fqp,iaT(»  m s l o * 5 * , 3 } 

ACCEPT  31.HSLUW 
T ypE  1800 

1800  FORMAT  ( * hF  AST  * 

ACCEPT  31,MFAST 
Type  1Sj3 

1900  FORMAT ( » TMAXs  '.51 

ACCEPT  JI.thax 
PNNTcn  »0.0 
PRMT  C2) sTMAX 
PRM r (31 =HFAST 
phmt(4) =0.001 
E*E*1  ,0/F|_OaT(N1 
00  2000  1*1, N 

MOOT ( I) sEWE 

wcn«wic(n 
200a  Continue 

n s w = i 

call  rkgscprmt.w. woot.n,  ihlf.fw.o*,  auxi 

N S W s N S W “ ! 

DU  2200  1*1, NF 
no  2100  Jal.NS 
AL(I, J) =0.0 
2103  CONTINUE 

2200  CONTINUE 

DO  2 430  1*1, NS 
00  2300  J*1,MF 
H(I,J)=0.0 
2300  CONTINUE 

240J  CONTINUE 

00  6230  ITER*1, NITER 
00  2830  1*1, NF 
00  2730  J * 1 , NF 

osci, J) *o(i,  j) 

IFCI-J1  2600.2500,2600 
2503  01(1,11*1.0 

GO  TO  2700 
2633  OI(I,J)*0.0 

2703  CONTINUE 

2803  CONTINUE 

CALL  LEaTlF(OS,NF,NF,  13,0I,3,'^KA«EA,  IER) 

00  3230  I*J ,NS 

00  3130  J*1,NS 

S U M K * 0 , 0 

Ou  3000  Ksl.NF 

3UML  * 3 . 0 
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2903 

300U 

3103 

3200 


3303 

3404 

3503 

3600 


3703 

3303 

3^03 

<1003 

4133 


4200 

4300 

4400 

4503 


4603 
47  00 


4300 

4900 


SUmL3SUML+0I (K,L) «CQ(L,J) 

CONTINUE 

SUMK*SUMK*ecI,K) * SUML 
CONTINUE 

A ( I , J) sA  ( I , J5 -SUNK 

CONTINUE 

CONTINUE 

00  3603  1st  1 NF 

00  3500  J*1,MS 

3UMK  =0 , 0 

00  3400  Kal,NF 

3UM|„  = 2 • 3 

Ou  3330  L»l, NS 

SUMLaSUML+CO (K ,LI *A (L, J) 

cuntinue 

SUMK=SUM«+ni f I , K 3 *SUNl 

CUNTINUE 

CN(I, JTsSUNK 

continue 

CoNT  InUE 

00  4003  T * 1 , NF 

00  3900  Jal.NF 

SUM*\  = 3 , 3 

00  3330  K»l,NF 

SUML  *3  • 3 

Og  3700  L = 1 , N S 

SUML aSUML ♦CO(K, LI *R (L, J) 

CONTINUE 

sunk 3 SUM* +0 1 (I  ,M  *SUML 
CONTINUE 

0(1, J) aOCl. JlfSUMK 
CONTINUE 
CONTINUE 
00  4400  Ial,NF 
DO  4300  Jal.NS 

SljH*0  1 0 

00  4230  K a 1 , Np 
SUMaSUM+OI(I,K)*COCK.  J1 
continue 

ALfI* J) a^LCI. JI+SUM 

CONTINUE 

CONTINUE 

00  4700  I a i , NS 

00  4630  J a 1 , NS 

*8(I,J)»A(IfJ) 

CONTINUE 
CUNTINUE 
00  4900  I a 1 , Np 
00  4820  Jal.NF 
OS  (I , J] *0 (I , J) 

CONTINUE 

CONTINUE 

CALL  £IGRFCA3,NS.7.0,UA,Z, 1 ,WKAREA, IER1 
CALL  EIG»F(OS,NF,  13,0, UO,Z,  t,'4K area,  t£9) 
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00  5102  1 = 1,  ms 

PRINT  5002,  (A (T, J) , J=l,NS) 

503J  F0RMATC'  i.Gl2.5.9nx,Gi2.5n 

5 1 3 J CONTINUE 

PRINT  <52^0 

5 g u j F 0 R M a T f / / / / ) 

00  5430  T = 1,ns 
phInT  53*2.  UACn 

53*3  FORMAT  ( * '.G12.5.2X,  # ♦ J ' , 2 X , G 1 2 . 5 ) 

5423  CONTINUE 

PRINT  52*0 

IF(NF.ld)  55*0. 553*, 5720 
55*3  00  5622  1 = 1, N F 

PRINT  53*0.  (0(1, J) . J=1  ,NF) 

5o*0  CUMTINUE 

GO  TO  60*2 

5730  00  5633  1=1, NF 

PRINT  5020,  (n  ( I # J) , J=1  , t*) 

5633  CONTINUE 

PRINT  52'5,3 

00  5903  1=1, NF 

PRINT  5023,  (0(1, J) . J=ll ,NF) 

59*3  CONTINUE 

(,230  PRINT  52*0 

00  b 1 03  1 = 1 , N F 
PRINT  5330.00(1) 

61*3  CONTINUE; 

6113  0Q  6163  1=1, NF 

00  b l 33  J=l,NS 
C0(I, J3  s C N ( I , J ) 

6133  CONjlNyE 

6163  CONTINUE 

6^33  CONTINUE 

00  6603  1=1, NF 
00  6533  J=1,NF 
0S(I,J)=0(t,J) 

IF(I-J)  6430,6303,6430 
6320  OI(I,I)al.0 

GO  TO  6530 
6430  01 (I, J)«0.3 

6533  CONTINUE 

6633  CONTINUE 

CALL  LEQTlF (PS.NF.NF, 11,OI,0.WKAREA, IER) 

00  7903  ITER=l, NITER 

00  7030  1*1, NS 

00  6900  J=l,NF 

SUMK*0t0 

00  6803  X=1,NS 

SUML=0.0 

Ou  6702  L=i,NF 

SUML  = SUMI_  + 0 (K,L)*OI(L,J) 

6703  CONTINUE 

SUM* = SUM* ♦ A (I ,K) *SUML 
CONTINUE 


6330 
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BN (I , J) sSUMK 
69 30  CONTINUE 

7304)  CUNTINUE 

7 l p»(d  00  7400  1 = 1,  NS 

00  7303  Jsl.NF 
3um=3.3 

00  7203  Ksl,MF 
3UM=SUH+8 [I, K) *01 (K, J) 
7230  CONTINUE 

HCI,J)sH(I.J)i-SUM 
7330  CONrlN|j£ 

7430  CONTINUE 

7530  00  7703  1=1, NS 

DO  7603  J=1,NF 
B(I, J] *3N(I,J) 

7&30  CONTINUE 

7700  Continue 

7300  CONTINUE 

00  9203  1 = 1,  N 
1 1 * I “NS 
DO  9100  J = l , N 
JJsJ-NS 

IF  (II)  7903,7930,3500 
7900  IF ( J J)  3300,8033,8303 

8000  SUN  a 3 # 0 

00  8103  K s 1 , NF 
SUH3SUN  + H(I,K) *AL  (X, J) 
8130  CONTINUE 

IF(I-J)  8330,3233,8300 

8 23'J  TW  ( I , I ) a 1 ,3-SUN 

TR I Cl,  I) » 1.0 
GO  TO  9130 

3300  TR(I,J)s-SUM 

TR I Cl , J 3 >3.0 
GO  TO  9100 

8430  TW Cl , J) s-HCI, JJ) 

TRI (I, J) =H(I, JJ) 

GO  TO  9100 

8530  IF(JJ)  8600,8630,8700 

8630  TR(I,J)3*LfIT,.n 

TRI Cl, J) >-ALCII, J) 

GO  TO  9100 
8730  SUN30>9) 

00  8803  K 3 i , NS 
SUM*SUM4,AL(II,KJ*HfK,JJ) 
8330  CONTINUE 

IFCl-ji  9030,8903,9030 
8900  TRCI,I)*1,3 

TRI Cl, 13  >1 . 0-SUM 
GO  TO  9100 
9030  TH(I,J)33,3 

TRI (I, J) 3-SUM 
9130  CONTINUE 

9200  CONTINUE 
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Ou  9 4 ;j <j  1 = 1,  MS 

XI (£] a3.3 

OU  93355  J=l,M 

XI  CI)=XI  CD+TRfl,  J)  *WIC(J) 

93»5'.J  CQMy  INijE 

4 4(i  j CONTINUE 

00  963?  1 = 1,  S|F 

1 

EtA  (I) =3 .3 
00  9533  J=t,N 

ETA  CT)  =£TA  tn+T3(It,  J]  *WIC(J1 
953>1  CONTINUE 

P^Cty  CONTINUE 

ENS=l  <{,/FL0AT(N3) 

OU  97y0  1 = 1,  NS 

w U 0 T (I] a E M 5 
9733  CONTINUE 

PR'1T(l]=0.0 

PR*<T  (2) aTMAX 
Pt?MT  (35  aNSLQW 
PRMT (4]  =3,331 
N S S s 1 

CALL  RKGS (P9MT, XI, WoOT.NS, IMLF, FSLO# OSLO. AUX) 

NSS=MSS-l 

EwF* 1 ,3/FuOAT (NF) 

00  9303  1=1 ,NF 
MOOT (I) »EWF 
9aP5j  CONTlNUg 

PHMTCiJ"0.0 
PRMJ  (2]  aTMAX 
PRMT C 3 5 aHFAST 
PRMT (41=0,^31 
NSF=l 

Call  RKGSCpRUT,ETA.wOaT,NF, I«LP, FFST . OFST,  AUX] 

NSF=N3F-1 
CNST=1, 3/416. 67 
TYPE  9900 

9933  FORMAT!'  GIVE  9AU0  RATE  ',S) 

ACCEPT  11.I8AU0 
130J0  TYPE  13103 

13130  FORMAT!'  WHICH  STATE  WILL  9E  PLOTTED?  *',S) 

ACCEPT  ll.IPLT 

IF  CIPLT-21]  13140,13113,10110 
10110  IR= 1 1+2* (IPLT-21) 

AMI Nil ,0630 

AMAx=_AMlN 

00  10130  I«i, NSW 

PUT (I) a0,0 

00  13120  J = 1 , 3 

JP4a  J + 4 

JP7aJ>7 

PLTCnaPLT(I]4FA(IP,JP4)*wS(JP4,  D+FACIR,  JP7]«WS(JP7,  II 
13120  CONTINUE 

PLT(I] =-CNST*CpLT(I) +FA (IP, 17] *WS(17, I ] *F A (IP , 1 9 ) * W S ( l 9 , 


1 


w 
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13130 

10140 


10300 


1 0300 
1 


i j5aa 

13600 


10  7 00 
10830 


13930 


1 1000 
11100 


11330 

11330 

U«30 

11420 


1 I) ) 

IF  CaMIN.GT .PLT  ( II 3 amin  = PlTCI) 

IFCPLTtn .GT.AMA*)  AMAX=PLT(I) 

CONTINUE 
Gq  TO  1 1«30 
D(J  11100  1 = 1. N5F 
WH(I3  =0.0 
TIM£sTF5T(I) 

DO  10300  JTal ,MSS 
JF  = JT 

IF  CAe$CTlME-T3L0CJT)  J-t  .0E-35)  10330. 13300, 10333 

CONTINUE 

GO  TO  ia<J00 

TINTig.a 

GO  TO  10800 

NMiaMSS.l 

DO  10700  JTal.NMt 

jf=jt 

IF (TIME-TSLO  t JTn  10700,10700,10500 
IF (TIHE-TSLOfJTtlJ ) 10600,10700,10703 

TNUMaTIME-TSLOIJT) 

ToENOH«TSLOCJT*n  -TSUOf  JT) 

TINT”  TNijM/TOENOM 
GO  TO  10800 
CONTINUE 
DO  10900  J=1,NS 

WH(I)»WHfIJ+TBI (IPLT. J1*CTINT«(XISCJ. JF+IJ-XISCJ. JP) )♦ 

ixiacj.jEn 

CONTINUE 

DO  11300  J=t,ME 

JJaJ+NS 

WH(I)  swHfD+THICIPLT#  JJ)  *£T43  (J,  I) 

CONTINUE 
CONTINUE 
AMlN=i .0E30 

AMAXa-AMIN 
DO  11300  1=1, NSW 

IF  (AMIN.  GT.  W3C!  PLT,  in  AM  ZNaUS  C I PLT , I) 

IFCWSCIPLT.  I)  .GT’.AMAXJ  AMAXaWSflPLT,  I) 

CONylNyE 

DO  11300  1=1, NSF 
IFCaMIN.GT.WH(I) ) AMIN=WH(I) 

IE  (WH  (t)  ,GT,  AMAX3  AMAX=WH(I) 

continue 

Ou  11400  1=1, NSW 
PLTCIJ aWSCIPLT. n 
CONTINUE 

call  inittcibauo) 
call  binitt 
call  dlimy camin, amax) 
call  nptscnswj 

CALL  xE«M(U) 

CALL  r FPM(fl) 

CALL  CH£CK(TW,PLT) 


CALL  QSPLAV  (TW.PLT) 

CALL  FRaM£ 

IF  (IPLt-21)  U440,  t 1^651, 11  at,  3 
1 1 440  CALL  LINE (34) 

CALL  QLIMY  (AMIN,  AMAX3 
CALL  NPTS(MSF) 

Call  check  (TF3T,whi 
CaLL  C°L0’’  t TFST  , hhi 
11460  CALL  ANMOOE 

ACCEPT  l#IANS 
TYPE  lisag 

l 1 5j0  FORMAT  ( • MORE  PLCT5?  ',S) 
l 1 62C1  ACCEPT  l.IAMS 

IF(IANS-'Y')  1 1 7 51  id  , 10000, 1 1 7 03 
11730  I F ( I AMS-  ' N ' ) ltaP'3,ll9«0,lt8fl'J 
llao'd  TYPE  40,lAN$ 

Gu  to  u fiao 
li^uo  type  120^2 

12 000  format (•  rum  another  case?  '.si 
12100  ACCEPT  l, IANS 

IFCiaNS-'Y'I  12200,100,12200 
12200  IF  ( I A N S « ' N ' ) 12300,12400,1230  0 
12300  Type  40, Ians 
Go  TO  12100 
12400  STOP 
END 

i *******su0routime  f’r  page  2 **«*«**« 

SUBROUTINE  FW (TIME, r.loot) 

DIMENSION  ,WPOT(N) , F A (20,20) 

CQMmON/bLKI /f», n 
OU  200  I=1,N 
HOOT(I) = 0.0 

00  102  J a 1 , N 

WOOTCI)  aWQOTCn+FAfl.  J5  *W(J) 

100  CONTINUE 

200  CONTINUE 

return 

ENO 

1****#**SUSR0UTIN£  FSLO  PAGE  3******* 
SUBROUTINE  FSLOCTIME, XI.NQOTl 
DIMENSION  XI (NS) ,WDQT (NS) , A (7, 7) 
CQMM0N/5LK2/A, MS 
DC  200  1=1, MS 
HOOT  (I) =0,0 
00  10a  J = 1 , NS 

HOOT  (I ) sROOT(n+A(t , J)  *XI  (J) 

100  CQNtINijE 

200  CONTINUE 

RETURN 
End 

1 *******SURRnuTINE  FFST  PAGE  4******** 

subroutine  ffst(time,eta,*dot) 

DIMENSION  ETA (NFl , MOOT  CNF) ,0(13,13) 
CQMMON/aLK3/0, MF 
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00  2*u  1*1, NF 

MOOT  (I)  = 0.0 
OQ  1 00  J*1.NF 

W0OTCIJaW0OT(n+O(I,  J)*£TaCJ) 

100  continue 

203  CONTINUE 

RETURN 
END 

1 ****r***SUBRnUTIME  Ow  RAGE  5******** 

SUBROUTINE  OWCTIME.w.wOOT, IHLF.N.PRMT) 

DIMENSION  W(N)  , WOQT  (N)  , PRMT  C5)  , WS  (20. 5311  , TW  (53H 

CONNON/OLK^/WS, TW, NSW 

OU  100  1=1. N 
WS(I,NSwJ C 1 3 

103  CONTINUE 

Tw  CnSwJ *Tihe 

IF (NSW.GE.50t)  PPMT(5)*1.0 

N S W * N S W * 1 

RETURN 

END 

l ******»*SUB ROUTINE  OSLO  PAGE  «,******** 
subroutine  oslo (time, xi , woot, ihle, n, prmd 
DIMENSION  XI  CMl ,W0OT(N) , PRHT ( 5 1 , X I S ( 7 , 50 1 ) , T$LO (50 1 ) 
COMMON /0LK5/X I S, TSLO.NSS 

DO  100  I 3 1 , N 
XlS(IfNSSj«Xl(IJ 
100  CONTINUE 

TSLO  CNSS) -TIME 

IF  (NSS.0E.53n  PPMT(S)*1.0 

MSSsNSS+1 

RETURN 

ENO 

1**  *SU0PQUTINE  OFST  page  7******** 

Subroutine  ofstctime.eta.wdqt, Ihlf.n.prmtj 

DIMENSION  ETA(N) ,hOOT(N] ,PRMT(5),ETAS(l3,50in 
1TFST (5ui) 

COMMON /dKb /ETAS,  TEST,  NSF 

DU  103  1*1, N 

ETAS (I ,NSF) =ETA  (1} 

100  continue 

TFSTCnSF) sTJME 

IF (NSE.GE.50l ) PRMT (5) * 1 .0 

NSFaNSF+1 

return 

END 
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• **«****t,*ppog9am  v a s m s p page  t ********* 

DIMENSION  A(7,7),xlC(7),XEQ(7).P(71,ASe7,7),0xlC(7lt*K(7] 
DIMENSION  41 (2.2) ,*2(2,5) ,*3(5,2) ,*4(5.5) , * 4S (5 , 5 ) , A 3 I 
1(5, 5) , A4C (5, 5), nxslC(2),DXFIC(5),A3(2.2),X(7),PRMT(5), 
2QERX(7),AUX(9.7),X8AR(2),ZaAR(5),ZwOL(5),XS(7,501),TXC 
3T(521),X6APS(2,531),Z8A9S(5.50l),T3AR(50i),Z'ir,LscS#501), 
4fwGL(5an  .PLT(5Pn  ,CTM?(3)  ,EFDI^T(53n  ,neLI^T(53n  , 
5XWGLSC2.5«U)  ,XC  (2,5*1)  , TxC(S0t)  # C C 1 8 ) 

COMPLEX  U(7) , Z(1 . 1) , VT, SG, Il,US(2) ,UF(S) 

PEAL  K A , ke, *f 

CQmmON/3LK3/KA  , KE, X F 

COmMQN/0LK9/H,  0 

CUHH0N/8LKI/ASAT,3$aT 

CQMM0N/3LK2/C 

COM'iOM/BLK3/VREF  , PIN 

cuf,M0N/aLK4/xs.  TxcT^'sxcr 

COMMQN/0LK5/ X3 APS , Z6 APS, TgAR,  NSXPP 

CUMM0N/3LKfe/ZwGLS , TwGL .nSZWGL 

COMrlON/BLK7/ZBAR 

COMMON /&LK 1 3/XC, TXC,NSXC 

EXTERNAL  FXCT,FXflR.FZwGL.OXCT,OXPR,OZ*GL,FZwGLC 
DATA  TOQP/5.0/, TOQP/3.5/ 

Data  xO/i ,2/,xo/i .*/. xp/0.25/ 
data  TA/0.G6/,TE/o'.5/,TF/1 ,0/ 

120  type  2512 

200  popmatc'  changing  any  parameters?  ' , si 

320  ACCEPT  4*0. IANS 

430  FORMAT  C a 1 3 

IF (IANS-'Y  ;)  503,300,500 
530  IF(IANS-'N;)  600,6100, S30 

620  TYPE  7 0 0 , I * N S 

730  FORMAT  ( ' ??',A1,'??  TYPE  AGAIN  *,S) 

GO  10  300 
830  TYPE  900 

630  FORMAT  ( ' WHICH  0N£?  ',$) 

ACCEPT  1000, IANS 
1303  FORMAT  (Afl) 

IF (IANS-'H')  1100,2500,1120 
1100  IF(IAhS-'O')  1200,2903,1230 

1203  IF (IANS-;TDOP;)  1300,3000,1303 

1303  IF  ( lANS-'TQGP')  1400,3203,1403 

1400  IF  ( I A N S - * X 0 ' ) 1500.3430,1500 

1500  IF ( I ANS-'XO')  1600,3600,1600 

1600  IF (IANS-'XP')  1700,3800,1700 

1700  IFClANS-'KA;)  1900,4000,1900 

1803  IF (IANS.;KE;5  1900,4200,1900 

1903  IF(IANS-'KF')  2000.4400,2003 

2330  IF  (IANS- 'T  A **)  210  0,4600,2100 

2100  IF (IANS-'TE')  2203,4800,2203 

2203  IFCIANS-'TFM  2233.5300,2230 

2233  IF (IANS-  'AS  AT' ) 2260 , 520  0 , 226  3 

2260  IF  (IANS-;3SAT;)  2330,5400,2300 

2303  TYPE  2433, IANS 

2400  FQRmaT(»  HEY  DUMMY  ',*4,'  is  NOT  ONE  OF  ThE  CHOICES'/'  TRY 
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1 AGAIN') 

GO  TO  830 

25c?a  TyP£  2633. H 

2633  FORMAT!'  OLO  Ms  ',513.5/'  N£W  Hs  '.Si 
ACCEPT  2733, H 

2733  FORMAT  !G1 

GO  T°  5633 

233^  type  2900,0 

2933  FQRmaTC'  OLO  0=  ',513.5/'  NEW  0=  ',51 

ACCEPT  2703,0 
GO  TO  5603 

3033  TYPE  3 t 0 3 , T Q 0 P 

3133  FORMAT!'  010  TDQPs  '.G12.5/'  NEW  TOQP*  ',5) 
ACCEPT  2703, TOOP 
GO  TO  5o32 

3230  TyPc  3333, TQOP 

3330  FORMAT!'  OLO  TGGPs  ;,G12.5/'  NEW  TOQP*  ',5) 

ACCEPT  2703,  TGiOP 
GO  TO  5630 

3430  TyPE  3500, XQ 

3530  FORMAT!'  OLO  XO*  ;,G12.5/'  NEW  X0  = ',S1 

ACCEPT  2733, YO 
GQ  TO  5630 

3630  TyPE  3703, XU 

3730  FCPMAT!'  OLD  X0  = '.G12.5/'  NEW  XG=  ',51 

ACCEPT  2703, XQ 
GO  TO  5630 

3330  TyPE  3 9 a 0 , X P 

3R3.J  FORMAT!'  OLO  XP=  '.Gta'.S/'  NEW  XP  = ',51 

ACCEPT  2703, XP 
GO  TO  5630 

40  30  TyPE  4103, KA 

4130  FORMAT!'  OLO  X A t ' , G 1 2 . 5 / ' NEW  K A = ',S1 

ACCEPT  2700, KA 
GQ  TO  5630 

4230  TYPE  4303, KE 

4330  FORMAT!'  OLO  KE  = '.G12.5/'  NEW  KE*  ',51 

ACCEPT  2700, KE 
GO  TO  5630 

440a  TYPE  4503, KF 

4500  FORMAT!'  OLO  KFs  '.Gta'.S/*'  NEW  K F = '‘,31 
ACCEPT  2700, KF 
GO  TO  5633 

4630  TYPE  4700, TA 

4700  FORMAT!'  OLO  T A ■ '.G12.5/'  NEW  TAs  ',$) 
ACCEPT  2703, TA 
GO  TO  5630 

4300  TyPE  4900, TE 

4930  FORMAT!'  OLD  TE=  '.G12.5/'  NEW  TEs  »,|) 

ACCEPT  2700, TE 
GO  TO  5600 

5030  TyPE  5100, TF 

5130  FORMAT!'  OLO  TF  s ''.G12.5/''  NEW  TF*  ' , SI 

ACCEPT  2703 , TF 
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GO  TQ  5o*3 

52*J  TyP£  5330, A5 AT 

53*3  FURMaTC'  OLD  ASATs  ',G12.5/'  NER  ASATs  •,$) 
ACCEPT  270U,ASAT 
GO  TO  56*3 

54*0  TyP£  550*,B5AT 

553d  EQSNaTT'  OLD  d S A T = ' , G 1 2 , 5 / ' ME*  8SAT*  ',51 

ACCEPT  a7<jafB5AT 
5o*J  TYP£  5723 

57*M  E 0 P M A T C * CHANGING  anymore?  ;,S) 

53*0  ACCEPT  4*0. IANS 

IE  Cl  A NS  - 'Y')  5900,300,590* 

593d  IE  (IAMS-'N')  o*4j*  , 6 1 *0 . 6*3* 

60*0  TYPE  7*3. IANS 

GO  7 0 56*3 

6i0j  type  0230 

62*3  FORMAT  ('  ***FnTER  LOADING  INFORMATION***'/'  VT* 

ACCEPT  6333.VT 
633(3  FORMAT  T2G) 

TYPE  642* 

64*3  FORMAT  ( * SG=  ',*) 

ACCEPT  630*, SG 
ILsCCJNJG  (SG/vn 
PlLsPEALf ID 
AI  IL*  AIM  AG  T ID 

RVT3REalTVT) 

AIVT=AIMAGTVTJ 
AnuM=xQ«AIIL-RVT 
OENQNsXQ*RTL*AIVT 
TH  = A7AN2  f ANUH, DENQH) 

ALPHA  = aTAN2  (AIVT.Rvn 
6E7ASATAN2  f 4riL,RILJ 

VMAGsCABS  (VT1 
AlLMAGsCABS  (ID 
ANGsAL°H A-TH 
VOrVuAG*COS  (ANG) 

VQSVMAG*SIN(ANG) 

ANGsBETa-TH 
AI0«AILMAG»C0S (ANG1 
AIQ»AILMAG*SIN  f ANG) 

X IC  C 1 D =V0+XP*AIQ 

x!C(5) ■xicci)+(xn-xp) *aio 

XICC2)=KF*XIC(5)/TF 
XIC (3) =TH 
XIC  (a) * i .0 

XlCC6)sC<E^SE(XlC(5)))*XlCC5) 

XICC7) *V0-XP*AIQ 
VREFsxIC(6)/KA+VMAG 
PINaREALfSG) 
type  650* 

65*3  FQRMATC'  GIVE  FINAL  AQhITT A N C E ' / ' Y2*  ',3) 

ACCEPT  272*, Y2 
TYPE  660* 

66*3  FQPMATC'  V2s  *,S) 
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ACCEPT  2733, va 
type  67a3 

6733  POPMATC'  IS  THERE  A CHANGE  IN  VREF?  '.$1 

6833  ACCEPT  400, IANS 

IF  ClANS-'Y')  6903,7130,69.30 

6900  IF(IANS-'N')  7030,7300,7000 

7300  TYPE  733.  IANS 

GO  T°  6303 
7100  TYPE  7200 

7203  FORMAT  ( * NER  VPEF=  ',$) 

ACCEPT  27a3,VREF 
7300  DO  7930  1=1,6 

xEOCDsxicm 
7900  CqNtIN|j£ 

C(D*i.a/Toop 

CC2)S1.3*C*0-XP)*Y2 
C(3)=(kO-XP)*Y?*V2 
C(«) *1 .0/TF 
C (5) =*F*C  (9) 

C (6] 80.5/h 
C(7) =Y2*V2 
C (3] =1 .0/TE 
C C9) =l ,0/TA 

C(10)  = (l.Z-Xt>*Y2W1.0-XP*Y21 

C(U)  = (1.2-XP*Y21*C(7)*xP 
CC12)=XP*XP*CC7]*C(7) 

C(13) =1 ,0/TQOP 
C(l9)  = (xO-tP5  *CC7) 

CCl5)=l.n+CXO-XP)*Y2 
CClb)*Ctl«)/CCl51 
c (17) =0 . 5*C  ( l 6) 

CC18J  =KE+C<A*KF/TF] 

TYPE  7910 

7410  FORMAT  ( * 00  LINEAR  ANALYSIS?  *.S) 

7420  ACCEPT  4O0.IANS 

IF(IANS-'Y;J  7430,7450,7930 
7432  IFflANS-'N'J  7940,14450,7940 

7440  TYPE  700. IANS 

GO  TO  7423 

7450  DO  3200  I TER  = l , 53 

I CON V s 1 
DO  7600  1=1,7 
00  7530  J = 1 , 7 
A(I»J)*0«0 
ASCI, J) =0.0 
7500  CONTINUE 

7600  CONTINUE 

VT8AR sSOR T CC  C10)  * CXEQ  Cn*XEQ(n+XEQ  (7)  *XEQC7))  *2.0*0(11) 

t*CXEQC7) *C0SCXE0C3))-XEQC1) *SINCXEQC3) ) )*CC12) ) 
P(1)sC(1)*CCC2)*XEQC1)*CC3)*SINCXEQ(3))-XE0C5)) 

PC2) =CC4) *CXEQC2)-C(5) *XEQC5) ) 

P(3)«377.0*£-XEQ(4)*1.0) 

P(4)=C(6)*C-PIN/X£Q(4)*C£7)*CXE3(l)*CQS(XEQ(3))*XEQ(7)* 
ISIn (XEQ (3) ) ) *0* (XEQ (4) -1 .3) ) 
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pc5)*c?ai*fCl<e+SEC^e,3(5ji)*xEQf5i-xenc<»n 
P(6)=C(9)*CKA*i:-XEQ(2)+C(5)*YEQ(5)+VTa4R-vR£f)*XEQ(6)) 
PC7)aCCl3)*C-CCl<»)*COSfxESf3n+C(l5]*XEaf7)  J 
A(l, lJ*-C(t J *C(2) 

A (1 ,3) *-C (1 ) *C(3) *COS  (XEQ(3) ) 

A(1,5)*CC1) 

*(2,2) *-C(a) 

* (3,5) *C (a) *c  (5) 

* (3, a) =377.3 

A(a,t)=_CCfe)*Ct7)*C0.S(XET(3)) 

AC«,3)sCC6)*C(7)*CXE'3Cl)*3lNCXEQf3n-XE3(7)*CQS(XEn(3))) 

A(a,a)s-ccfe)*(PlN/(x?'j(4)*xET(/oi>0) 

A(a,7}s-C(b)*C(7)*SiNCXEQC33) 

*(5,5) *-C(8) *(KE+(R3*T*xGQ(5)fl .0) *SE (XEO (5) ) ) 

* ( 5 , 6 ) =C ( 8 ) 

A (6,  l)s-(C(<9)*KA/VTBAH)*(CC10)*XEQCl)-C(ll)*SIN(X£Q(3))) 
A (6, 2) =C (9) *K A 

A(b,3)  = (C(9)*KA/VT8A9)*(C(ll)*(XEQ(7)*SI,‘J(XeQ(3))+X£Q 

1(1) *C03 (XEQ  (3) ) ) ) 

A (8,5)  s-C(9) *K  A #C (5) 

A (6,6) =-C (9) 

A (6,7)  s-(C(9)  *KA/VT8AR)  *rC(l?)  *XEQ(7)  +CC1  l)  *COS(XE'3(3)  ) ) 

A(7,3)«-Cn3)*CCl41*SlM(X£Q(3n 

A(7,7)s-C(13)*C(15) 

00  7820  1=1,7 
OQ  7700  J s l , 7 
A3(I,J)=A(I,J) 

7700  CONTINUE 

7800  CONTINUE 

Call  LEQTlECA,  1 ,7, 7,P,0,WK,  IER) 

00  8130  1*1,7  . 

IF (ABS (P (I) ) -2.0001 ) 8000, 8000, 7 P00 
7900  ICONVsg) 

8000  XEQ(I)sXEQCI)+P(I) 

8100  CONTINUE 

IF (ICO MV)  8200,6200,3400 
3203  CONTINUE 

type  3300 

3300  FORMAT ('  £0  POINT  CONVERGENCE  FAILURE') 

STOP 

8400  PRINT  8500 

8500  FQRMATC'  a * * OLO  EQUILIBRIUM**#') 

00  8700  1*1,7 
PRINT  86^0, I,XTC(I) 

8600  FQ9MaT('3',2X,'XEQ(',I1,')=  ;,GI2.5) 

3703  CONTINUE 

PRINT  8800 

8800  format c * ***new  equilibrium***'’) 

00  6900  1*1,7 
PRINT  3600,  I.xEQCI) 


8903  CONTINUE 

PRINT  9300 
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PRINT  9400. (AaCl.JI , J=l,7) 

940y  FORMAT  ( ' '.7(2X,G12.5) 1 
9500  CONTINUE 

OU  9600  I=t,7 
OXICfnsXIC(I)-XEQ(I) 

9600  CONTINUE 

PRINT  q j?!Z 

9700  FOrmaT  ( ' ***HKJE#.R  INITIAL  CONDITIONS***'*] 
00  9900  1*1,7 
PRINT  9500, 1 , 0 X I C (I) 

9600  FORMAT C '0* . 2* . 'OXIC  (' , 1 1 , ']  = '. 012.5] 

9900  CONTINUE 

00  9960  1=1,7 
no  9930  J=l,7 
4(i,j)SAScr,j) 

9930  CONtInijE 

9960  CONTINUE 

CALL  EIGRF(AS.7,7,0,U,Z,t,WK.IER) 

PRINT  10000 

13000  FORMAT  C'  ***SYSTEM  EIGENVALUES***'] 

00  13200  1=1,7 
PRINT  13130.UCT] 

10130  FQRmaT  f ' '.2X.G12.5,'  +J  ',G12'.5] 

10230  continue 

DO  10400  1=1,2 
00  10303  J = l . 2 
A1 (I, J] =A(I, J] 

10303  CONTINUE 
10430  CONTINUE 

00  10600  1=1.2 
00  10500  J*l,5 
JJaJ+2 

*3(1, J) =A(I,JJ] 

10500  CONTINUE 
10600  CONTINUE 

00  10300  1*1,5 
00  10700  J = 1 . 2 
11=1*2 

A3CI,J]=ACII, J] 

13700  CONTINUE 

10800  CONTINUE 

DO  11000  1=1,5 
00  1090O  J*  1 , 5 

II»I+2 

JJaJ+2 

A4 ( I , J] s A ( I I , J J) 

*4S(I| J)*Aa(i, J) 

10900  CONTINUE 

11000  CUNTINUE 

00  11400  1=1.5 
00  11330  J ■ 1 » 5 
IF(I-J)  l 1 200 ,11100,11200 

11100  A4 1 (I,  I] =1.0 
GO  TO  11300 
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liars'  A4i(i,J)»a.a 
11300  CONTINUE 
H4a?  CONTINUE 

CALL  LEQT1PCA4.5.5.5.  Aai.vj.wK,  !£»} 

DO  11603  1=1.2 
DO  11733  J = 1 » 2 
SUHlag.a 
DU  11600  <=1.5 
3u*2=y .0 
DO  U50O  L = 1.5 
SUM2s3UH2*42fI,U  *A4l  (L.X) 
l 1500  CONTINUE 

SUMi=3UHl+SU^2*A3C<, J) 

116^0  CONTINUE 

acci.  j)  *n  (i. 

11700  CONTINUE 
11800  CqnTINUE 

DO  12200  1=1.5 

00  12100  J=1 .5 

SU.^iaa.O 

DO  12000  <=1.2 

SUM2=0,3 

DU  11900  - = 1.5 

SUm2  = S'JM?  + A4I(I,1)*a3CL,K) 

11Q00  continue 

SUMl=SUHl+SUM2*A2  C< , J) 

1200O  CONTINUE 

AaC  CI,J) =A4S  a . J) +3UN1 
12100  CONTINUE 
12200  CUNTINUE 

00  12300  1=1.2 
DXSICCn«OXICCI) 

12300  continue 

00  12600  1=1,5 
11=1+2 
SuM 1 . 3 
Da  12500  <=1.2 

S u M 2 = 3 • 3 

00  12900  L = 1 . 5 

SUM2  = 3UM2  + A4I  CI,U  *A3CL,<) 

12400  CONTINUE 

SUMi  sSUMUSUM2*DXlC  (K) 

12500  CONTINUE 

OXPIC(I)  =OXIC(in+SUMt 
12600  CONTINUE 

PRINT  12700 

12700  FORMAT ( ' ***SLO'*  SUBSYSTEM  MATRIX***') 

00  12503  1=1,2 

PRINT  R400 . (AQCI.J) . J=l,2) 

12800  CONTINUE 

CALL  ETGRF  CA0,2,2,0,US,Z, 1 ,WK , IER) 

PRINT  12900 

12900  FQRmAT  ( ' »**SLOW  SUBSYSTEM  EIGENVALUES***'’) 
00  13000  1=1,2 
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PRInT  1 0 1 0(* , US  ( I ) 

133lD  CONTINUE 

P H I N T 1 3 K1 0 

1 3 1 FORMATf*  «**SLOW  INITIAL  CONDITIONS***') 

OU  133510  1 = 1.  a 

PRINT  i3?53f IfOXSIC(I) 

13200  FoaMATf'f^.a^.'OXSICC'.Il.'la  012.5) 

13300  CONTINUE 

print  ijaaa 

13430  format  C'  ***uncorhfcted  fast  subsystem  matrix***') 

00  13500  1=1.5 
PRINT  94?<J.  (A4S(If 
1 35*j0  CONTINUE 

CALL  EIGRF{A4S.5.S.0,uF#Z, 1.WX.IER) 

PRINT  1360O 

13600  FORMAT ( ' »**UNCOPRECTEO  FaST  EIGENVALUES***') 

DO  13700  1=1,5 
PRINT  i3l00fUFCI) 

13700  CONTINUE 

PRINT  13300 

13800  FORMAT ( * ***CORRECT£0  FaST  SUBSYSTEM  MATRIX***') 

DO  13900  1=1.5 

PRINT  <?a00.  £A«C£I#J)  , J=i,5) 

I39w0  CONTINUE 

CALL  EIGRFf A4C. 5.5.0, UF, Z, l,WK, IER) 

PRInT  1400O 

14000  FORMAT  f • ***CQR«ECT£D  FaST  EIGENVALUES***') 

00  14130  1=1.5 
PRINT  10100, UFfl) 

14133  CONTINUE 

PRINT  14200 

14230  FQRMAT(»  ***FaST  INITIAL  CONDITIONS***') 

DO  14400  1=1,5 

PRINT  14330, I .OXFIC(I) 

14300  FORMAT C'3',2X, ;0XFIC C ' , 1 1 ,') = ;,Dl2.5) 

14430  CONTINUE 

14453  Type  ia502 

t45<J0  FORMAT  ('  ***PREPARING  FQR  SIMULATION***'/'  FAST  STEP 

1 5IZE=  '.*) 

ACCEPT  2T30,HFAST 
TYPE  14510 

14513  FORMAT ( ' SLOW  STEP  SIZE*  ',S) 

ACCEPT  2703 , HSLOW 
TYPE  14520 

14520  FORMAT  C ' Tmay*  ',$) 

ACCEPT  2700, TMAX 
PRMTC4) =0.301 

Type  14530 

14533  FORMAT ('  IS  THIS  A FAULT?  ' , S) 

14540  ACCEPT  400. IANS 

IFfIANS-'Y;)  14550.14570,14550 
14553  IF(IANS-'N')  14560,14750,14560 
14563  TYPE  700. IANS 
GO  TO  14543 
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1457?  7r  = E \4530 

1453/T  FO».‘:AT(*  * * *Ph  GCOOY-A  FAULT***'/'  GIVE  V?  DURING 

lFA'JLT  ' . 31 
ACCEPT  270?, Y2F 
TYPE  14590 

1459?  FC9MATf'  GIVE  FAULT  OIjPaTION  ',5) 

ACCEPT  2700, TCLP 
CTMPCilsCCP) 

CTNP(21 «CC3) 

CT*PC3)sC(7) 

Ct^P  (41 =C(1?1 

c TMp (51 sC  ( 1 1 1 

CTMP(b)sCC12l 

CTMp(7)*C(lUl 

CT*P(a)=C(l51 

DO  i4fc?3  1=1/7 

KCDaXICtn 

1 4 6 vj  P C i j N t I N u E 

C ( 2 } = 1 ,0*  ( X Q - X P ) *Y2F 
C (31 =3.0 
C (7) =0.0 

C(lOls(1.0-xP*Y2F)*(l.?-XP*Y2F1 

cci  i)»a.9« 
c da)  =«?.a 

C ( 1 4 ) 30 , 0 

C (151 =i .0+  (xq-xpi *y2f 
p««t ([) *0. a 

PRMT (21 =TCLR 
PRMT  (31 s P F A 5 T 
Erf=l  .0/7 .3 
OU  14650  1=1,7 
DEPX (II sEw 
1465?  CONTINUE 
NS  X C T = 1 

Call  RKGS(PRMT.XfOERX,7,IhLF,FXCT,OxCT,AUXl 
CC2) =CTMP(1) 

C (3) *CTMP  (23 
C C 73  s C 7 N P (3J 
C (101 =CtmP (41 
C (1 11 sCTMP (51 
C U 21 «CTMp(bl 
C(141sCTPP(71 
C ( l 51 sC  TMP (31 
00  14700  1=1.7 
XIC(I) »X(I5 
TYPE  *,X(I) 

14700  CONTINUE 

PAUSE  'IC' 

GO  JO  14350 

14750  00  14300  1=1,7 

X (I) =XIC(I1 
148^0  CuntInuE 
14350  PRMT(11=0.0 

PRMT (2) > 7 M A X 
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PK“>T  (3)  =HFAST 
Ew =1.0/7. 3 
00  14903  1=1,7 
D£PX  (I) SEH 
14R'30  CONTINUE 
N s X C T = t 

CALL  RKGS(PHMT,X,0ESX,7# IHLF.FXC7J0XCT, AUXi 

NSXCTsmsXCT-I 

PkmT(1)=0.0 

PRMT  (2) sTmax 

p«mt  (3) aMSLQW 

OQ  15303  1=1,2 

D£RX  Cl) =0.5 

x9ak  cn  sxiccn 

15000  CONTINUE 

0 u 15353  1=1,5 
11*1  + 2 

ZBArCI)  =XIC(H) 

15053  CONTINUE 
NSXbRs  i 

CALL  RKGSCPRMT,X8AP,DERX,2,IH(.F,FX8R,ax3R,AuX) 
MSX6PsNSX8R-i 
00  15133  1=1,2 
X8AMfI) sXIC  (I) 

15130  CONTINUE 

Oq  15153  1=1,5 
11=1+2 

zbar(I)  sxiccin 

15150  CONTINUE 

CALL  50LVFA CXSAP, Z8 AR) 

PRMT  f l ) =0.0 
PRMT (2) =TMAX 
PRMT  C3) =HFAST 
Do  15200  1=1,5 
11*1+2 
OERX (I) =0.2 

zwglcD  =xiccin-zsARci) 

15200  CONTINUE 

NSZrGl=1 

CALL  RXG3  CPRMT,ZWGL,0ERx,5, IHLF.FzwGL.OZWGL, AUX) 
NSZwGL=N3ZWGL-1 
T yP£  15300 

15330  FORMAT**  ***PREPARING  FqR  PLOTS***'/*  GIVE  BAUD  RATE  ', 
IS) 

ACCEPT  15400, I8AU0 
15400  FQRMaTCI) 

TYPE  15410 

15410  FORMAT  C * PLOT  ZERO  OROER  APPROX?  *,S) 

15420  ACCEPT  400. IANS 

IFCIANS-'Y*)  15430.15450,15430 
15430  IF(IANS-'N')  154a0, 15460, 15443 
15440  TYPE  732. IANS 
GO  TO  15420 

15450  Call  SNGPL7CX5,TXCT#NSXCT,X8ARS,Z0ARS,T0AR,NSX5R,ZWGLS, 
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lTwGL#NSZrfGl.l*AU0.7,2,5,NSXB».T8A!») 
t 5 4 a J TYPt  15535, 

15530  FORMAT (•  00  CGPRECTIONS?  * , 5) 

15600  ACCENT  400. IANS 

IK(IANS-#Y*3  15703,  15900,  1 5 7 c» 0 
157  33  IF  ClANS-'N*)  15803.  19400,  153^3 
158.J51  TrP£  7 3 j , IANS 
GU  TO  15800 
1 5 Q J 3 E F D I N T (1 ) =3.3 

no  16333  I=2,NS2WGL 
in  1 = 1-1 

EFrHNT{  I)  sEFOINT  (Ikin  +5J.5*  (TWGL  (T)-TWGl  C T M t ) ) * CZ*GLS  (3. 
lp+ZHGL5(3,If-M)) 

10000  CONTINUE 

EFUSUMsEFO  IMT  CMSZWGL) 

OEL  IMT  C 1 ) =3.3 
Nax^^l = N S X H R - 1 
OU  1 7 b 3 3 I = 2.NSZWGL 
I M 1 = I _ l 

DO  17500  J = l,2 
JT  = I- J+  1 
TIMEsTMGLCJTI 
no  162510  jXst.wgxBR 
JFIN0=JX 

IF  ( A0S  ( T ImE-ToAR  f J X ) ) «•  l , 3 E - 3 5 ) 16300.16330,16203 

1 b 2 J 3 CONTINUE 

Go  70  16600 

16300  IF ( J - i ) 16400,16400,16500 

lo«00  ziBRKs^B ARS (l , JFINO) 

GO  TO  175O0 

16500  Zlt3HKl=Z5Ac’S(l,JFINOl 
GO  TO  17500 

1 66 J0  DO  16900  JXsl.NSXBMl 
JFIN0=JX 

IF(TIME-T0AR(JX))  1693^,16900,16700 
16700  IF (TIME-TBAR CJX*l) ) 16830,16900,16900 

16300  TNUhsTIME-TQARfJX) 

TqENQMsTBAP  CJXM1-TBAR  ( J X ) 

TINt  = tN(jH/TOENOM 
GO  TO  17200 
1 b 9 0 0 CONTINUE 

OlFF*TIME-T8ARf JFINO) 

TYPE  17003, TINE, TbARlJFlNO) ,OIFF 

17300  FORMAT!*  ',3f2X,E15.8n 

pause  'corect* 

TYPE  17130 

17100  FORMAT ( * INTERPOLATION  FAILURE  OURInG  CORRECTION*) 

STOP 

17203  ZlNT*TlNT*(ZBARSIl,jFlNO  + n-  ZBARS(l,JFINn))+ZBARS(l,JFlN 
10) 

IFCJ-i)  17333,17300,17400 
17300  zl3rRKszlNT 
GO  TO  17503 
17430  Z1*RKi»ZINT 
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1 7 500  CUMTINUE 

Fk»SIn(218RK)  * f-l  ,51  + COS  CZwGLS  Ct  , I J 1 J +COS  ( Z l 3°*)  * S I N j Z WG 
1LSC1»  I)  ) 

FKlaSXNCZlBRKi) * (- 1 .3  + CQS (ZwGLS C 1 , IMl J1 ] + COS C Z 1 SR K 1 ) * S I 
lNCZWGLSCl.  1*1)3 

DEL  I NT  Cli  = DEL  I NT ( IMl) +3.5* (TWGL (I) -TtfGt (IHI ) ) * (FK  + FK1) 
17600  CQNTInUE 

OELsUNsOELINT (NSZwPL) 

CNSTUCCO  *CC3) 

CNST2»C («]  *C  (5) 

AKls-CNSTl tOELSUH+C ( 1) *EF0SUM 
X8AR  (15  sXIC  CD  +AK  1 
AK2aCNST2*EFOSUM 
X0AR (23 sXIC  (2) *AK2 

Type  », x*ar ci 3 , xpa«  (2) 

TYPE  17  75*43 

1770C*  FORMAT  ( • PRINT  CORRECTION  RESULTS?  ',53 
177<j3  ACCEPT  4530,  IANS 

IF(IANS-'Y')  177516,17712,17753b 
1772b  IF(IANS-'N')  17709,  1 7790  , 177539 
17749  TYPE  70S , IANS 
GO  TO  17703 

17712  PRINT  17720, EFOSUM.OELSUM 

1 7720  FQRMATC'  '.G15‘.a,5X,Gl5.83 

OU  17760  1*1,  NSZWGL 

PRINT  1 77453  , TWCL  (13  , ZWGlS  (1  , T3  , Z^GLS  (3.  1 3 , EFO  I N T ( I ) , 
IQE|_INT(I) 

17740  FURNATC'  '.G15.3,4(5X,G15.3)  3 

1 7 7653  CONTINUE 

00  17780  1*1,  NSXPR 

PRINT  17740(T8ARCI3,ZBARS(1,T) 

17780  CONTINUE 

17790  pause  'coric; 

Oq  17800  I=1,NSZWGL 

XWGLSCli I) *-AKl-CNSTl*CELINT(I3+C(l) *£FQINT(I) 

XWGLS (2, I) =-AK2+CNST2*EFDINT(I3 
17800  CONTINUE 

Oq  17900  1*1,5 
Z6ARCI3  sZBARSCI,  13 
17900  CONTINUE 

PRMT (l)a0.0 
PRMT(2)*TMAX 
PRMT (3) sHSLOw 
Oq  18000  1*1,2 
OERX (I) *0.5 
18000  CONTINUE 

NSXbR  * i 

CALL  RKGS(»RMT,XBAR,0ERX,2, IMLF,FxBR,0X8R, AUX) 

NSX8R=NSX8R-1 

NSX8M1*NSXBR-1 

00  19100  I=l,NSZWGL 

TIM£=TI»GLCI) 

OQ  18100  JT*l,NSX6P 
JFINQbJT 
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IFCABS(TINE-T3AR(JT))-i,0E-35)  1*2^3,  192'1J,lBl03 
tdljM  C QN  T I NUE 

GO  TO  18433 

1 d 2 .j  0 r>Q  133^3  J = l . 2 

XC  ( J , n sVwGlS  ( J,  n +X9a®S  C J,  JF  I NO  3 
1 a 3 0 3 CONTINUE 

gu  rn  19100 

13433  oa  i <9 7 V* :j  JT=  ( , NSXtj*  l 
J F I '•  0 = J T 

iFCTiME-TaAPfjTn  1 37,10,  is7j0,  tasoo 

la*)?  TP  c TIhE.T8AR  f JT+1  ) 1 136113,1(3730,18733 

166J0  TNUMsTIME-TdAR  ( JT) 

T0EN0MsT9A®  CJT+n  - To AH  (JT) 

T IMT  = TMU'm/  TOENOM 
GO  TO  lgGa? 

1 6 7 j y continue 

OIFr  = riME-TdArt  CJF  I"JUT 
TYPE  1?O00,  TIME,  Ttf  Ai?r  JFJNO)  ,0  IFF 
PAUSE  'XC; 
type  lsa^o 

16*30  FUP-IATC'  INTFRPUL  ATXOn  FAILURE  at  xc') 

STOP 

loPwfl  Oq  193  33  J = 1 , 2 

YdW  = TlNT»(YdAFSCJ,JFlNn+lJ-XRAPS(J,JFINL)n+'XBAPSCJ,J 
IF  INC) 

XC  ( J , I ) a X ,( G L S ( J , I) +X®R 
19033  CONTINUE 
19130  Continue 

NSXC=NSZWGL 
CO  19*313  IS1  . N S Z ’w  G L 
T x C ( I ] sTWGL(T) 

19230  CONTINUE 

PkHT(1)!?,D 
PPM  T (2)  sTMAX 
PRM  T (3)  =HFAST 
Oq  19303  1=1.5 
OERXCI) *0.2 

Z«GLCIJsxICCI+3)-2ba9s(I,i) 

19330  CONTINUE 
NSZwGq* i 

Call  RKGS (PRMT, ZWGL.0E9X.5, IHLF.FZWGLC.OZwGl. AUX) 
NSZWGLaNSZ’PGL-l 

CALL  SNGPLT (XS, TXCT.NSXCT, XC. ZRARS, TBAR.NSXC.ZWGLS. TWGL, 
1NSZWGL, I a AUO. 7 ,2,5, NS  X 09, TxC) 

TYPE  19310 

19310  FORMAT  ( * WANT  ADDITIONAL  ITERATIONS?  ;,S) 

19320  ACCEPT  400, IANS 

IFCIANS-'Y*)  19330.15903,19333 
19333  IFClANS-'N')  19343,19403.19343 
19343  TYPE  700, I ans 
GO  TO  19320 
194q3  TyPE  19533 

1 95 J3  FORMAT ( * RUN  ANOTHER  CASE?  ; . $ ) 

19633  ACCEPT  400, IANS 
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IFUANS-'Y')  19700.  100,  1 7 59  4 

1 9 7 y 3 IF(I*NS-'M*)  19803,19900,19330 
TYP£  7«3.IA‘)S 
GO  TO  19631 
\99ya  STOP 
E NO 

l*#****j»*-*PL0CY  0 A T A PAGE  ?********* 

block  Data 

St»L  K A , *E  , k F 
CUMHON/^LKlyASAf.BSAT 

CO^MOM/aLKS/v A ,*E,KF 

CUMP.0N/dLK9/H,  n 

Oa^a  aSa^/o.00U23/,8SAT/0.30«1/ 
data  H/s.a/,n/2.e»/ 

haTa  KA/P5.;)/,KE/-0.i9«5/,KF/y’.  16/ 

EnO 

• ****»*»**c,ijMCTION  S F PaG£  3********* 

Function  sfcakgi 

CUMMnN/6LKl/ASAT,6SAT 
SEsaSAT  * E X p C^SAT  * A 9 G 1 
Pt^URM 
E,*0 

1 <***«*•  *<SL-nR0UTiM£  F x c t Page  a********* 

SuflPOuTP'E  FXCTCTIUE,  X.OERX) 

D I m e M S I G |k)  X (7)- , DERX  (7}  , C £ 18) 

REAL  KA.KE.^F 
COMMON /0LK3/C 
CU,nsON/rJLK3/VPEF,PIN 
cgmmon/bl* « /*  a , ke  , k f 

CU^f-iOM/bLKP/H,  0 

Cx3  = cascx  (3)1 
SX3*SINCX  £3n 

VLT  = SGRTCCfl0)*(X(t)*XCnTXC7]*X(7ns2.y*Cfll5*CX(7]*CY3 
l-X£i)»SX31*C(121) 

OERX  (i)*C  U ) * C-C  C2)  *X  (l)-C(3)  *SX3  + X £5)  1 

ne»X  £2)  «C  (a)  *(-x  (21  tC(5j  *x  (51 ) 

OtRX  C3) =377.1* ( x (ai -1  .1} 

f)£RX  £«3  sC  (6J  * (PIM/X  (<»J  - C (7)  * (X  (ll  *CX3+X  (71  *SX3l  -0*  (X  (4)  - 
1 l.Bl 1 

OtRX (51 ■C(8lB(-(KE*seCX£5lll*X(51+X£611 

HcPX  (61  sC  (9)  *(KA*  CX  (2)  -C  (5)  *X(51  -VLT*'/9EF)  -*  C 6 1 J 

OERX(7)aC(131*(C(14)*CX3-C(l5)*X(71) 

R£  TUR.'l 
Eno 

j .*  *******SUP«OUTINE  FX0R  Pa^E  5********* 

SUBROUTINE  FXBR (TIME, X0AR, OERX 1 
niMfc^'3IQN  xaAK(2l,OERX(2l,C(l81,ZBAR(5! 

COHMON/aL<2/C 
CQMMON/0LK7/79 AR 
CALL  SOLVFA(XBAR,ZPaR] 

nERX(ll»C(n*£-C£21»XBAR(ll-C(31*SlN(Z5AR£lll+ZB*R(311 
D£RX  (2l  »c  (4)  *C-XfUR(21+C(51  *ZBAR(31 3 
R£  TuRM 
E.'iO 
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■ **«r»r*#***SijNROUTlME  F Z N G |_  P 4 fj  E 6«**-»****** 

Su^PnoTr^E  FZhfiLIT  tme#z*gl,0ERx) 

ci  MENS  ION  ZWGUC5)  ,D£R X (5) ,C( 181 ,X8ARS (3,5*1 ) .Z^ARS (5,5*1 1 , 

1TBaR(5<)1)  »Z^AR(5)  »X8AR(2)  , Z ( 5 1 

B£At  K A , K £ , K F 

CO^-'Ort/ciLKp/C 

COMMON /OLK 3 /VRFF , PIN 

COMMON /bLK 5/ Xo  APS,  ZRARS,  T-3  AS  , -ISX3R 

C(jM,-tON/aLKF/K  A , AE  , KF 

COMMON/aLKR/n.O 

nu  i c»y  jTsi  ( MSYiit? 

J F I N 0 s . T 

IFCaBS(riMF-TbA«f  jT)  1-1  .*£-351  15*.  15.3,  1*® 
lo*  CONTINUE 

go  rn  33P 

15*  r,0  2«a  1 = 1.2 

XbAWtI)aX5A!?$CT,JFlNn) 

20*  CONTINUE 

G(j  28  3 1 = 1,5 
Z8AR(Il=ZbA«S(r,JFTNO) 

2S*  CONTINUE 

Go  Tn  75* 

3u*  Ns:<bNi=r,Sx8s-i 

no  450  JTsi  jMSXbvj, 

JFsJT 

IF  (TINE- T3 AH  f JT1  1 «50 , «5« , 35* 

350  IF  (TINE-T3AW  (JTfin  403  , 'AS*  , «5* 

400  TNUM*TIME-T6AH(JT) 

TOEhOfUTPAP  (JT*11-TBAR(JT) 

TlNTaTNUM/TcENOH 
GO  TO  bOO 
450  CONTINUE 

OIFFsTlNE-TfeSARf JF) 

TYPE  5*0. TIME, T6ARCJF) ,QIFF 

5uo  format  ( * *,3(2x,ei5.an 

pause  'fz*gl' 
type  55o 

53G  FORHAT(»  INTERPOLATION  F*IIU°E  IN  FZ^GL'l 

STOP 

60*  Dq  650  1 = 1,2 

xeAR(I)sTlNT*(X3AK8(I,JF-*-ll-Y3ARSCI,JFll*X8APS(I,JF) 

65*  CONTINUE 

Go  7P0  1=1,5 

Z8AH(I)sTlNT*(ZBARS(I,JF+ll-Z8ARS(T#JFl)+ZRAPSCI,jF) 

70*  CONTINUE 

75*  Oo  3*0  1=1,5 

Z(IjsZRAR(I)t-ZWGL(I) 

300  CONTINUE 

czwi=coscz(in 

SZR1=SIM(Z(H  1 

VuTaSCRT  (C  ( 1 * ) * ( X S A R III  *<3*R(11+Z  (5)  *Z  (51  1 ♦2«0*C(H)  * ( Z (5 ) 
1 «CZ*1-XBAR (l )»3ZW1) *C(12)) 

D£RX ( 11  =377.0*ZWGL  (2) 

OERX  (2)  «C  (6)  *(PIM/Z  (2)  -C(7)  *(XRA9(11  *CZxU2£5)  *S7’Rl  )-0* 


W 
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1 ZWGL (2)  3 

OE*X  (3)aC(8)*(-(KE+gE(2(3)))*Z(3)+Z(fl)3 

0£RX  C4)aCCR)*CKA*CX8AR(2)-C(5) * Z ( 3 ) * VREF- VL T) -Z ( a ) ) 

DERXC53  aC(13)  * CC (1 «) *CZw l -C ( 1 53 *Z(5) 3 

return 

END 

» ******«**surroutine  o*ct  page  7********** 

SUBROUTINE  QXCTCTIME.X.DERX,  IHLF.NDIN.PRMT) 

OIM£NSlQN  X(7)  ,0ERX(7)  , PRmT (5) , XS (7 , 501 ) ,TXCT(501) 
COMHON/BL^a/XS. TXCT.NSXCT 

no  iee  i » t . 7 

XS(I,NSACT)«X(I) 

U0  CONTINUE 

txCt  (NSXCT) *T Im£ 

IFCNSXCT.GE.5ai)  PRMT(5)si.0 
NSXCTaNSXCT+l 

return 

EnO 

1 *******»*su8routime  oxaR  page  a********** 

SUBROUTINE  OX8B (time. xBaR,DERX. IHUF.no IM.PRMT) 

0 1 MENS  IQN  X8AR(25 ,0ERX(2) ,PRMT(5) ,X8ARS(2,501) #TBAR£50l) , 
1ZBAR(5) ,ZPaRS(5.5CT1) 

CONMON/BUK5/X0ARS, ZBARS, T8AR.NSX8R 

CUMM0N/eUK7/Z8AR 

CALL  SOLVFACXBAR.ZBAR) 

DO  100  1=1.2 
X8ARScr,NSX8R)=XBAP(I) 


180  CONTINUE 

Du  200  1=1.5 
ZBARS(I,NSXBR) «Z8AR(I) 

280  CONTINUE 

T0AR  (NSX8R) =TIM£ 

IF(NSX8R.GE.501)  PRH T (53=1.0 

NSXBRsNSXBRtl 

RETURN 

EnO 

i *********SUBROUTINE  OZwGl  page  r****««**» 

SUBROUTINE  OZRGL (T IME, ZWGL.OERX.IHLF.no IM.PRMT) 

DIMENSION  Zwr,L(5)  ,0£RX(5)  ,PRMT(5)  ,ZWGLS(5,50l)  ,TWGL£S0l) 
C0MM0N/BLK6/ZWGLS.TWGL.NSZWGL 
DO  100  1=1,5 
ZWGLS(I.NSZWGL) *ZWGL(I) 

100  CONTINUE 

TWGL(NSZWGL)*TIME 

IE  (NSZwGL.GE.50l)  PRMT(5)«1.0 

NSZWGL«NSZWGL+t 

return 

EnO 

i»#«******3ubroutine  solvfa  page  i0**«*«***»* 

SUBROUTINE  SOLVFA (XBAR, ZBAR) 

DIMENSION  X8AR(2) ,ZBAR(5) ,C(18) 

real  ka.ke.kf 

C0MH0N/8LK2/C 
CQMMQN/BLK3/ VREF. PIN 


. 


W 
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CUH.-tON/aL*  1/ASAT.BSaT 
CQMMON/FILXa/KA.KE^F 

CUMM0N/5LK9/H, D 

0ITaZ3ARCl) 

00  i 05d  IT£R  = 1,50 

ANUNs-PIM  + C (7J  * (XbAB (11 *COS  CO  I T J *C (17) *S IN (2. 3*0  IT) ) 
0eNU»'«cc7j*C*B*R(n  »StM(0lT)-C(l6J  *C03  (2. 0*0  IT)  ) 
OELsanum/OENOm 
niTaOIT>OEL 

IF (ABS(oEL) -0.00011  3S0,300,  i net 
130  CONTINUE 

TYPE  20'J 

2 jn  FuRMAT(#  CONVERGENCE  FAILURE  ON  OELT  A-50LVF  a * ) 

STOP 

30p  zaAam*oiT 

Zb  AH  (5)  aC  (lb)  *COS(ZdAR(l) ) 

Z6AR(2)ai  .c? 

VTBARaSQRTfCf 10) *(XBAR(l) *XBAR  f 1) +ZBAR (5) «ZBAR (5) ) +2.0*C  C 
11  1)  * (ZB  ARCS]  *CUS(Z8AR(1)  )-XRAR(lJ  ASlNfZBARCD)  ) + C (12)  ) 
0IT*zaARf3) 

00  400  I TER  * 1 , 50 

ANUMa-(C(l8) ♦SE(OIT) ) *0IT*KA*  CX8AR(21 +VREF-VT3AR1 
0ENUMaCC13)+faSAT*0lT+l ,0) 

OELa ANUM/OENOM 
0 1 TaO I T *OEL 

IF (AB$ (DEL) -0.00*1 ) 800,800,400 
430  CONTINUE 

TYPE  500 

5j0  FO»NAT(#  CONVERGENCE  FAILURE  ON  EFO-SOLVFa') 

STOP 

8 JO  ZoA«(3}*0IT 

ZBAk (4) «* A* (X8AR(2)-C(5) *ZBaR (3) +VBEF-VT8aR) 

return 

ENO 

1 ****»**#**SU9RUUTINE  fz^glc  PAGE  1 !••******•• 

subroutine  fzhglcctime.zwgl.oerx) 

0 1 M£NS ION  Zi-r,Lf5)  ,OERx  (5)  ,C(ld)  »XBARS(?,50l)  ,Z8ARS(5,501 
1)#T8ARC501)  ,XC(2,50l)  ,TXCt50n  ,Z3ARI51  ,X(2J  ,Z£5) 

REAL  ka,ke.kf 

COPMON/8LK2/C 

COMMON /6UK 3 /VREF, PIN 

C0MN0N/8LK5/X8ARS,Z8AWS, T8AR, NSX8R 

CUMMCN/6LK8/KA,KE,*<F 

CUHM0N/6LXR/W,0 

CONMON/eLKl0/XC,TXC,NSXC 

00  100  JTal.NSXC 

JFIND«JT 

IF  CA9S(TINE-TXC(JT) )-l  .OE-05)  200,200, 100 
100  CONTINUE 

GO  TO  400 


200  00  300  1*1,2 

X (I) «XC (I.JFINO) 
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400 

530 

630 

703 

800 

900 

1300 

1100 

1300 

1300 

1400 

1503 

1633 

1 700 
1 80y 

1903 

2000 

2103 

2 200 
2300 

2403 


NS*CM1 sNSxC-1 

00  700  jTsl,MS*CMl 

JFINOsJT 

IF  (TIM£-TXC(JT5 ) 703,703,500 
IF(TlME-TXCCJT+m  600,700,703 
TNUMsTIME-TXC  (JT) 

T0EN0M»7^C  (JT+U-T7C  (JT) 
tINt*tnum/tuenom 
GO  TO  1000 
CONTINUE 

OlFF=TIME-TxC (JFINO) 

TYPE  600. TIME, TXC CJFINO) ,0IFF 
EU9HATC*  ',3(2*,E15.8) ) 

FaUSE  'FZRGLC' 

TYPE  900 

FQPMAT  C * interpolation  failure  IN  FZWGLC-XC*') 

STOP 

Ou  1100  1=1,2 

X(I)sTINT#(XC£I,  JFINO  + IJ-XCCI,  JFlNOn+XCtl,  JFINO) 
CONTINUE 

00  1300  JT= 1 , NSX8N 
JFINOsJT 

IF(a8$(TIME-T8AR(JT))-1,0E-O5)  1400,1400,1300 

CONTINUE 

GO  TO  16®0 

00  1500  1=1,5 

ZaAR(I)sZ8ARS£I, JFINO) 

CONjINyE 
GO  TO  23®0 
NSXbM13NSX8R-l 

00  1800  JT=1,MSX8M1 
JFINOsJT 

IF CTIME-T8AR(JT))  1900,1900,1703 
IF(TIME-T8ARfJT+l))  1800,1903, 1900 
TNUMsTIME.TBAR(JT) 

TqENOMsTBAP  (JT+l)“T0AR (JT) 

TINT*TNUm/TOENOM 
GO  TO  2100 
CONTINUE 

OlFF*TlME-TaAR (JFINO) 

TYPE  800, TIME, T0ARCJFINQ) ,DIFF 
PAUSE  'FZWGLC' 

TYPE  2®S® 

FORMAT  ( * INTERPOLATION  FAILURE  IN  FZWGLC-Z8AR') 

STOP 

Ou  2200  1=1,5 

ZSAR(1)=TINT*CZ8ARS(I,  JFIMO^-n-ZBAPSCI,  JFINO))  A-ZBARS 

1 (I, JFINO) 

CONTINUE 

00  2400  1=1,5 
ZCI)=ZBARCI)*ZWGL(I) 

CONTINUE 

CZW1«C0S(Z(11) 

SZ«l«SIN(Z(l)) 
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VLT«5QBT(C(l«)*CX(l)  •X(l)tZC5)»Z(5J  J*a.a*CUl)*C2(5J* 
lCZ"*t-X(l)  *SZw  1 ) +C ( 1 2) ) 

OE3x(ns377.0,ZWRLf2l 

DE**  C2)  «C(6)  * (Piw/z  (2)  -C(7)  * (x  (l)  *CZWl+Z  (5)  aSZM)  - 
10*ZWGL (2) ) 

0ePX(3)«C(8)*C-(KE*sECZC3ni*ZC3)+ZC«)  1 

DE^x  (4)  =C(<>)  * (KA«  (-Y  (2)-C  (5)  *Z  (3)  *VR£F-VLT)-Z  (4) ) 

OEPX  (5)  aC(131  « (C(l«)  *CZ-4l-C(151  *Z  (5J  V 

RETURN 

END 

J **********$U8POUTINE  SNGPLT  PAGE  12********** 

SUBROUTINE  SNGPLT (ES,TX#NSTpP,X8ARS#Z8ARSfT8AR,NSTPX8,ZWS, 

1 T*  , N5TP7.W  , XDAUO , N , NS , NF , NS TPZ8 , TxC ) 

01  MENS  I cr1  ES(N^4STPE),Tx(NSrPEl,XBARS(MS,NSTPX8]  ,Z8A»S(NF, 
IN5TPZH)  . Z'*4  3 (NF.NSTpZK)  , T 3 AR  ( NS  7»  Z« ) , TW  ( MsTPZ  W 1 , PL  T ( 53  11  , 
2Z4aT(5GJ) «rxC(M5TPx8) 

4yfl  FORMAT (At) 

DO  640fl>  I PUT* 1 , NS 
AMIMsl .3E3D 
AMAXa-1 ,<*£30 
00  6003  JTai.NSTPE 

IFCAMIN.GT.ES  CIPLT.J Til  AMINsESCIPUT. JT) 

IF (ESCIPLT, JT)  ,GT. AMAXl  AM  A X *ES ( I PL T , JT ) 

6JP3  continue 

Ofj  6100  JTsl.NSTPXB 

IF (aMIN.GT.X8APS(IPlT, JT))  AM  I N = X 8 APS  ( I PLT , JT ) 
IF(xaARS(IpLT# JT1 .GT.AMAX)  AMAXaX8 4RS ( IPLT, JT) 

61C10  CONTINUE 

Du  62510  JT«1,NSTPE 
P|_T(JT)=ES(IPLT,  JT) 

62cm  CONTINUE 

CALL  INITT(I8AU0) 

CALL  BINITT 

CALL  OLImYCAHIN, AMAX) 

CALL  NPTS(MSTPE) 

CALL  XFRM(4) 

CALL  Y F R m ( 4 ) 

CALL  CHECK  (TX  , PLT) 

CALL  OSPLAY  (TX.PLT) 

CALL  FRAME 

00  ©3SG  JTai.NSTPXS 

PLT(JT)8X3ARSCIPLT,JT) 

6303  CONTINUE 

CALL  LlNE(3tt) 

CALL  OLIMY (AMIN, AMAX) 

CALL  NPTSCNSTPX8) 

CALL  CHECK  (TXC. PLT) 

CALL  CPLOT(TXC.PLT) 

CALL  ANMOOE 
ACCEPT  4P0 , I anS 
6400  CONTINUE 

00  7600  IZHATel.NF 

ZMATf l)aZ0ARS(IZHAT, l)*ZHS(IZHAT,  1) 

00  6000  ITZwa2,NSTPZW 
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00  6500  ITXa*2.NSTPZe 
NTX8*ITX8 

IFCTWCiTZVO  .LE'.TBAR(ITX8)  ) GO  TO  6700 
6500  CONTINUE 

IF  ( (TW(lTZH)-TBAR(NTXB) )-l .OE-04)  6703.6700,6550 
6550  TYPE  6600 

6600  FORMAT  ( * ***SEARCH  FAILURE  0U9ING  INTERPOLATION****) 

STOP 

6700  ZlNTs  C ( TW(TTZM)-T8AH(NTx8«| ) ) * (Z8AR$  ( IZHAT, nTXB) -Z0APS 

t C IZHAT, NTX8-1) ) / (TBAR(nTx8)-TBAR (NTXB- 1) )) *ZB ARS ( IZHAT, 
2MTXP-I) 

ZHAT(ITZ^) =ZWSf IZHAT, ITZWj  *ZINT 
6800  CONTINUE 

IZhAt2*IZHAT+NS 
AMINsl .0E30 
AMAxs-l .0E33 
00  6900  JT»l,NSTPE 

IF (aMIN.GT.ES CIZHAT2, JT))  AM IN*ES ( I Z* AT2 , JT ) 
IF(ES(I*HAT2, JT) .GT.AMAX) AHAX*ES(IZHAT2, JT) 

6900  CONTINUE 

00  7000  JTM.NSTPZW 
IF(AMIN,GT.ZHAT(JT) ) AMIN.ZHATC JT) 

IP(ZHAT(JT1  .GT’.ahax)  AMAXaZHATI JT) 

7C00  CONTINUE 

Ou  7100  JT  a 1 , MSTPE 
PLT(JT)*ES(IZHAT2, JT) 

7100  CONTINUE 

CALL  INITTCI8AUQ) 

CALL  8INITT 

CALL  OLlNYfAMIN,AMAX) 

CALL  NPTS(NSTPE) 

CALL  xFRM(4) 

CALL  YFRM(a) 

CALL  CHECK(TX.PLT) 

CALL  DSPLAYCTX.PLT) 

CALL  FRAME 
CALL  LINE(3u) 

CALL  OLI^y (AMIN, AMAX) 

CALL  NRTS(NSTPZW) 

CALL  CHECK (TW , ZHAT) 

CALL  CPLOT(TW.ZHAT) 

CALL  ANMOOE 


ACCEPT  aOe.IANS 
AMIN.1 .0E30 
*MAx*-1.0E30 
00  7200  JT»1 .NSTPZ8 

IF(AHIN,GT.Z9ARS(IZhAT, JT))  AM IN*Z9 ARS ( IZHAT . JT) 
IF  (ZB  ARS  ( IZH  AT.  JT)*.  GT.AMAX)  AM AX *ZB ARS C IZHA T , JT) 
7200  CONTINUE 

Dq  7300  JT* 1 , N3TPZW 

IF(AMIN,GT.ZWS(IZHAT, JT) ) AMIN*ZWS (IZMAT, JT) 

IF  (ZWS(  IZHAT,  JT)  ’.GT.AMAX)  AMAX*ZWS(IZHAT,JT) 


PLT(JT)  »73AF?S(IZHA^.  JT) 
7«*c»vj  continue 

CAUL  INITTdBAUCn 
CALL  SINITT 
CALL  OLImY f AMIN, Ah  AX) 
CALL  NPTS  (MSTPZtjl 
•CALL  xp5MC«) 

CALL  YPH'-»(4) 

call  cnecKcreup.PLn 

CALL  OSPLAYCTdAP.PLT) 

call  fpame 

no  75^  jw.nstpz* 

P|_T(JT)  sZ-ASUZMAT,  JT) 
75*0  CONTINUE 

Call  LlNE(3a) 

CALL  uLIMYCAMI.n,  ANA*) 
CALL  M P T 3 C N j>  T P Z <0 
Call  CHECK  (T.-«,PLT) 

CALL  C°lOTCT^,pl1') 

call  am-icoe 

ACCEPT  4*4, IANS 
To'-'o  CONTINUE 

S t T j 3 M 

END 
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1*********PR0GRAM  mmsp  PAGE  t********** 

DIMENSION  H(3)  ,0EL(3)  ,XQ(3),VO(3) ,VQ(3) , XIC(20) , XP(3)  , 
1X0(3) ,TF(3) ,VR£F(3) ,PIN(3) ,Y(3,3),TH(3,3) ,PRMT(5)  , 

2X0QT  (20) ,X (20) ,AUX (8,20) ,XS (23,501 ) ,PLT (501) ,TS (501) 
3,0(3) , IA(3) , TE  (3 ) ,T00P(3) ,T00P(3) , VL  T (3) , X3A«(7)  , 

AP( 13) ,X9ARS (7,501 ) ,Z«ARS(t 3,501 ). TRAP (501) ,XC (7, 
5501),TXCC501),ZHGLS(13,501)  , 7WGL(50t)  . Z'JGL(13)  , 0ELPW(3)  , 
feAlNTl  (501)  , AINTP(5t,l)  , AIMT3(501)  , AIMTa(53l)  ,XWGLS(7, 
7501) , A (N 15(501) ,0PUP(3) ,ZOQT (13,501) , TOOT  (501) 

COMPLEX  TTC, YF(3,3) . YPO (3,3) ,SG(3) . VT(3) 

PEAL  IMaG,IO(3),IQ(3),KF(3),KE(3),KA(3),IDB(3,501),IQB( 
13,501) , IOw (3,501 ) ,IQW (3. 5«1) 

COMMON /bLK 1 /XO , XO, XP 
CQMMON/fjLKP/TQOP.TOOP 
COMMON /GLK3/T A , TF.TF 
CUMMON/yLKa/H.n 

C0MMCni/aLK5/ASAT,bSAT 
CQMriON/aLK6/»<A,KE,KF 
COMMON /bLK  7 /Y , TH 
COMMON /BLK 8 /BH , SH 
C 0 M H 0 N / B L K P / V P E F , P T M 
COMMON/3LK10/XS,TS,NXCT 
COMMON /bLK 1 1/VLT 
COMMON/iiLK  1 2/ZPAP 

COMMON /GLK 13/XPAPS,ZBdRS,T5AP,NXR« 

COMMON/ftLK ta/XC,  TXC  , NXC 

C0mn0;J/sL<15/ICPCT 

COMMON / bLK 1 6 /ZWGLS, TWGL.NZWGL 

CUMHON/bL^ 17/ZOOT, TOOT ,NOOT 

EXTERNAL  EXCT,OXCr,FXBAR,aXBAR,FZwGL,OZwGL 

DATA  (SG(t) , 1*1 *3) / (0.723,0.3703) , (1.63,0.0654) , (3.85, 

1-3. 1395)/ 

DATA  (VT (I) ,1*1,3) / (1.04,3.3) , (1.312,3.1648) , (1.022,0.0 
16292) / 

1 FORMA f (A  1) 

4 FoRMAT(Aa) 

n FORMAT  (II 

31  FqhMAT(G) 

32  F0RMAT(2G) 

SH»H(2)+H(3) 

8H=H (1) *Srt 

139  DO  200  1=1,3 

ITC«CONJG(SG(I)/vT(n  ) 

IMAG=CA8S(ITC) 

0ETA  = ATAM2(AIMAG(ITC) .REAL  (I  TC) ) 

VMAG»CA3S(VT(D) 

ALPHA  = ATAM2(AIMAG(VT(I) ) , REAL  (VT  ( I) ) ) 

ANUM*XQ(I)  »ArfiAG(ITC)-REAL(VT(r)  ) 
OENQMmXa(I)»REAL(ITC)+AIMAG(VT(I)) 

0£L ( D = ATAN2 ( ANUM, OENOm) 

ArlG*ALPH»-OEL  (I) 

VU(I) *VMAG*COS(ANG) 

VQ ( I ) «VMAG=StN(ANG) 

ANGsBETa-OEL(I) 


--  - ■ - 
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10 (I) =ImAG*COS(AMG) 

T HI)=IMAG*3IN(ANG1 

a ja  continue 

xiccn*i.5) 

DU  300  I 3 1.3 

xicci*ij«vQcn*xp(i}*io(r] 

120  CUNTIhUE 

OQ  <403  1 = 1,3 

xiccuni  *vgf  n+xoru  *iom 

aU0  CUNTlNUe 

OU  503  1=1,3 

XXC(I+4)«KF(I)*XICfI+ilJ/TF(I) 

520  CONTINUE 

Dq  1=1,3 

<iC(i  + ia)  ■(KE(n*sF(xic(iti  n n *x  ict  1*1  n 

63M  CuNTI^l.'E 

Ou  79*3  1 = 1,3 

xiC(iti7j«vo(n-xPfn*iaci) 

730  CONTINUE 

Ocj  3P3  1 = 1.3 

VHEFfI)a{XlCfI  + l«J /K A ( I ) ) +C A3S f VT  ( I ) } 

PIN(X)  aREALCSGTin 
8„0  CuNylMtjE 

XIC(8)a(HCU/8Hj#(fH(a)»0EL(2)*HC31*0eLC3))/9H-0EL(lJ  J 
XIC(9J=3.0 

XIC(13)=n£LC?3-OELf3) 
xicu  n *0.0 

TyP£  903 

<*00  FORMAT  {*  GIVE  STUDY  DURATION  #,Sl 

*CC£P  7 3 1 , Tma x 
TYPE  1030 

1303  FORMAT ( • GIVE  SLDw  sT£P  SIZE  ' , S) 

ACCEPT  J1,h$LQW 
TYPE  1130 

1103  FORMAT ( * GIVE  FAST  STEP  SIZE  '.$) 

ACCEPT  jl.HFAST 
TYPE  1220 

1203  FORMAT ( * IS  THIS  A FAULT? 

1323  ACCEPT  l.IANS 

IF(IANS-'Y»)  1400,1700,1400 
1403  IF (IAnS-*n*)  1500,2600,1500 

1502  TYPE  1600, Ians 

1603  FORMAT  C • ??*',Al,*'??  TYPE  AGAIN  *,*) 

GO  TO  1 3?0 
1703  TYPE  1830 

1603  FORM A T ( # GIVE  CLEARING  TIME  ;,S) 

ACCEPT  3 l , TCLE AR 
TyPE  1900 

1900  FORMAT  C * GIVE  NAM£  OF  FAULTED  V MATRIX  #',S) 

ACCEPT  4 , NN A h 
CALL  IFILE(20,NNAM1 
00  2030  1=1,3 

R£AO(20,321  (YF(I, J) , J«i,31 
CONTINUE 
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OQ  2333  1*1,3 
00  2203  J*I,3 
YCI,J)*CA85(YF(I,J)) 

Th(I, J) *ATAN2(AIMAG(YFCI, J) ) ,REAL(YF(I. J) ) ) 

IF  C I - J 3 2100. 22'03, 2100 
Y(J»X)«Y(I. J1 
TH(J,I)*Th(I,J) 

CONTINuE 

continue 

PRMT  C 1 ) 30.0 
PRHT(2) sTCLEAR 
PRMT (3) sHFAST 
PRMT (4) *0.001 
00  2400  1st, 20 

XOOT (II *0.05 
X(t3*xIC(I) 

CONTINUE 

NXCT=i 

CALL  RKG3CPRMT.X. XOOT,20, IHLF.FXCT.OYCT, AUY) 

NxCT3NxCT-l 

00  2500  1*1,20 

XIC(I) *X3(I,NXCT1 

TYPE  <*,XS(I,NXCT1 

continue 

PAUSE  »IC; 

PRMT  C j 1 *0 . 0 
PRMT (21 aTMAX 
PRMT  C 3 D sHFAST 
PRMT  (41 *3.001 
DO  2700  1*1,20 
XOOT(I) *3.05 
X (I)*XIC(I) 

CONTINUE 
TYPE  2800 

FORMATS  GIVE  FILE  NAME  OF  POST -DISTURBANCE  Y 
1 M ATR I X »,S) 

ACCEPT  4 , NNAM 
CALL  IFILE  (20 i NN AM) 

00  2900  1*1,3 

READ  (20, 32)  (YPO(I.J) , J*  l , 3) 

CONTINUE 

Oq  3200  1*1,3 

00  3100  J*  1 , 3 

Y(I, J)«CA8S(YP0(I, J}} 

TH(I, J) *ATAN2(AIMAG(YP0(I, J) ),REAl(Yp0(I,J))1 
IF  C I - J 3 3032,3100,3000 
Y (J#  I)  *Y  ( I » J) 

TH(J* I) *TH(I, J) 

CONTINUE 

CONTINUE 

NXCT*i 

CALL  RKGS(PRMT,X.XOOT,20, IHLF.FXCT.OXCT, AUX) 

NxCT*NXCT-l 

00  3300  1*1,7 
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xsArfcnsxiccn 
3 3'v0  CONTINUE 

Oq  3400  I*t , 13 
Z»JARCl)SXIC(Ii-7) 

3400  CONylN.jE 

PKMT  ( ^ ) 3*  .* 

PRMT (5) ST«4X 
PHH  T C3) =*SLOW 
PSMT  (4)  *0.031 
E*S=t  .3/7.0 
00  3533  I si,  7 
XQQT  C I ) sE.'S 
353 j CONTINUE 

NX*R3  1 

call  rkgs cprmt, xpa», xoot, 7 , ihlf,Fxrar,ox9ar, auxj 

NX0R  s^XgR - 1 
PRMrms3.3 
PHVf (21  3 T H A X 

PR  M T (31  a HF  AST 
PRMT (4}  33.331 
ErtF*l .3/13.3 
00  3600  1*1,13 
XOOT ( I ) «EwF 

Z*GL<U  ■*ICCI  + 7)-ZRaRS(I.  1) 

36i*0  CONTINUE 
NZWGLsi 
ICRCT=0 

CALL  RKGSCPRMT.ZWgl, XOOT. 13. IHLF.FZWGL.OZwGLiAUXI 

MZWgl=nZWGL-1 

TyPE  3700 

370J  FORMAT  ( ' GIVE  9 A U 0 RATE  #,S) 

ACCEPT  1 1 , 19  A'JO 
TyPE  3900 

3900  FORMAT!'  PLOT  ZERO  ORDER  RESULTS?  ',S) 

3900  ACCEPT  l , IANS 

IFtlANS-'Y')  4000,4203.4300 
aOO'O  IF  (I  ANS-'N  •)  4 100,4  303,4  1 30 

4100  TYPE  1630, IANS 

GO  TO  3903 

4203  CALL  GRaPHS(NXCT,NX3R,NX9R,NZwGL. IDAUO.XS,TS,XSARS, 

1T3AR,Z9ARS, TdAR, ZWGLS. TWGL) 

4303  TyPE  44Q0 

4403  FCRMAT(>  00  CORRECTIONS?  »,S) 

4500  ACCEPT  1 , IANS 

IF  (IANS- ' Y *)  4630,4803.4630 
4600  IF  (lANS.'N'J  4730,8700.4730 

4703  TYPE  1633, IANS 

GO  TO  4500 
4800  OELR  ( 1 1 *3 . 0 

OELRw  (1) s0.3 
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OELR  £2)»BH*Z3ARSf 1. IT)/H(l)+H£3) *Z8ARS£3.IT)/SH 
DELS (3) >Bh*Z8ARS(1 . ZT)/H(1)-M(3) *ZBARS(3. XT) /SH 
00  49a0  j=l,3 

ANGsTH  ' I , J) >OELR ( J) -OELR  (I) 

CANG=C0S(AMG1 
SANGsSIN  (ANQ) 

IOBCX,  IT) ■I08(I,rT)»VCI, J)*(ZBARS(J+ia,IT) *c ANG-X3 AR$  ( 
1 J+ l , IT) *SAN6) 

IU8(I#TT)«TQB(I,  IT)  *y  (I , J)  *(X8ARS(J+l,  in  ttCANG  + ZBARSC 
1 J*10,  IT)  >*SANG) 

4900  CONTINUE 
5000  continue 

00  5300  ITM.N7WGL 
IOw(I» IT) *0.0 
IQw  (I,  IT) *0.0 
TIMtsTWGL  (IT) 

CALL  INTPQL  (TIME, TBAR, MX 6R,JF, TINT) 

0ELC*t^nT*(zBAR5(1.  JF  + 1)-ZBA93(1,  JF) ) +Z8 ARS  ( 1 , JF) 

0£LO  = TlNTir  (ZflA9S(3.  JF+ 1)-ZBA0S(3,  JP)  ) +ZQARS  (3  , JF) 

OEL«  (2) *BH*06lC/M(1 ) + HC3) *OELD/SH 
OELR  (3) =BH*OELC/H(l3-H(2) *OELO/SH 
0£LRN(2)=8H*ZwGLS(t# IT) /H £1) *H  (3) *ZWGLS(3.  IT)/5H 
0ELRW(3) *BH*ZwGLS£1,  IT)/H(1)-H  (2)*ZWGLSC3, IT)/SH 

00  5200  J*  1 , 3 

5100  EQPbTIMT*£XBARS£J+1 , JF + i ) - X B AR S £ J+ 1 , JF) ) +xB ARS ( J* l , 

1 JF) 

EDPbTINT*£Z8ARS£J<M0#  JF*1)-ZBARS£J*10,  JF)  ) +Z3ARS  C 
lj*i0. JF)+ZwGLS£J+10*IT) 

EDPS=£OP-ZWGLS(J+i0, IT) 

ANG8=Th(I, J)+OELPCJ)-OElP(I) 

CANG*C0S(ANG8) 

SAMGs5IN(ANG8) 

ANGWsOELPN ( J) -DEL*w  ( I ) 

Tl*-i  .0*CUS  (ANG'rf) 

T2=sIN(ANGW) 

T3  = Tltt  .0 

I0W£I, IT) =IOW£I, IT)*Y£I| J)*£ (EOPS*CANG-enP*SANG) *Tl- 
1 CEOP*SaNG+EOP*CANG) *T2+ZWGLS£ J+10. IT) *C»NG*T3) 
IQW£l,XT)«IQW£If IT) +Y(I, J)*( (E0P*CANGtE0PS*3ANG) *T1* 

1 £EOP*CANG-EQP*SANG) *T2fZWGLSCJ+t3. IT) *SAnG»T3) 

5200  CONTINUE 

5300  CONTINUE 

5400  CONTINUE 

a in  t i (i)  *a'.0 

AINT2(i)*0'.0 

AINT3C1) *0.0 

AINT4  (1) *0.0 

DO  5300  IT*2.NZHGL 

IFCIT-2)  5500,5500.5600 

5500  FiMi*ZwGLS(2. 1) / (XBAPSC1# l)*XBARS£i  , 1) ) 

F2MUZWGLSC4, 1)  / CXB  ARS  £ 1 . I ) *XS  ARS  £ 1 , 1 ) ) 

GO  TO  5732 
F 1 M i *F  t 
F2M 1 *F2 


5600 
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573d  I M 1 = I T - l 

CALL  XNTPQL  CT*Sl  f IT j , T8A9, NX5R.JP , T I'm 

OMfiaTlNT*  CX8ARS  (t  , JF*l]  - X8ARS  Cl  , JM  3 + *3  (1  i JF) 

GMPS=aMR*QMS 

FUZWGLSC2,  IT)  /OMRS 

F2sZwGLS  C4, IT) /QMBS 

OT.TWGLCin-TrtGLf  INI) 

AINri(IT)aAINTtCl«l)*3.5*0T*fFtfFlNl) 

AINT2(IT) ■AINT3CIH1 } ♦J.5»OT#(F2*F2hI) 

AINT3CIT)  sAINT3(TM1  ) *3 . 5 *0T  * ( Z^GL  S ( 2 , I T J +Z^GLS  C 2 . I N 1)  ) 
AINT4  (IT)  ■ AIMT«  (TMl  ) +0.5*07*  f ZWGLS  C4,  IT)  + ZWQLSC4,  I'll  ) ) 
563^  CONTINUE 

AINT5Cl)s3.d 
CU  6403  ITs2,N7«GL 
INI *IT-1 

call  intpolctwgl  citj  , rn ar.nxrr.jp, tint) 

IFClT-2)  5903,5930 ,6t03 

59?d  FlHlsfcl.0 

00  6030  J * 1 , 3 

FlMiaFtMl+XaAR3CJ  + l,l)*lQW(J,n  + CZ5ARSCJ  + lR.l)+ZWGLSCJ 

1 + 13, i) ) *inwc j, n +zwgls (J+ia. n * tone J, n 

6o3j  CONTINUE 

GO  TO  62^3 
613  0 F I N j s F i 

6200  Fls0.3 

OQ  h303  J = l , 3 

£ U P s T I NT*  £ X 0 A R S (J  + l ,JF+l)-XBAR3CJ  + l,JFn+X5ARSCJ+l,JF) 
EOPsTlMT*CZ8APSCJ+l0#JF+l)-ZBA9S(J+l3#JF))+Z0A9SCJ+l3, 
i JF) +ZWGLS C J*1R, IT) 

AlOBaTlNT*fIOBCJ.JF+l)-l08CJ.JF))+I08CJ,JF) 

F 1 sF  1 +EQP  * I QW  C J , TT)+Eup*I0W(J,lT)+ZwGLS(Jt-l3.IT)*Al09 
6330  CONTINUE 

OTsTWGLCin-TwGLCIMn 
AINT5(IT)=AINT5(IMl)+0.5*OT*fFl+FiMl) 
fe400  CONTINUE 

Cl*  (0.5*Sh*PlN ( 1) / (H Cl ) *8H) ) * A INTI (NZWGL) 
C2=(3.5*PIN(2) /3H) * C-AlNTi (NZwGL)-H£3) *AINT2CNZ^GL) 
1/SH) 

C3s (0.5*pIN(3) 7 b H ) * ( - A T N T 1 CnZWGL) + h ( 2 ) * A I n T 2 ( N Z w G L ) 
l/SH) 

C4s3,5*Q (t) *SH*AINT3 (NZWGL) / CH(l) *RH) 

C5* (0.5*0 (2) /3H)*(-AlNT3 CNZWGL) -H (3) *AINT4CNZXGL)/ 
ism) 

C6S(0.5*O(3)/8H)*C-AINT3(NZWGL)-*'H(2)*AINT4(NZWGL)/ 

1SH) 

C7*-C0.5/8H)*AINT5(NZWGL) 

X8A* (I) aXIC (l) +C1+C2+C3+C4+C5+C6*C7 
00  6530  ITal.NZwGL 

Ti*“Ct  + C®«5*SH*PIN(l) / (H(l) * 5 H ) ) * A T N T 1 (IT) 
T2s-C2+(0.5*PIN(2)/8H)*(-AINT1(IT)-HC3)*AINT2CIT)/ 

1SH) 

T3a-C3+(3.5*pIN(3)/aH)*C-AlNTl(IT)*H(2)*AINT2(IT)/ 
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T5«-C5* (0.5*0  C2) /8W)*(-AINT3(IT)-m(3) *AINT4(IT) / 
1SH) 

To=-C6+C0.5*O(3) /3H) *C-AtNT3(IT)fH(2)  »AIMT4(IT)/ 
ISH) 

T7=-C7-(0.5/Bh)*AINT5(IT) 

XwGLS  ( l , IT)  = TH.T2  + T3<.r4.».T5  + T6  + T7 
6500  CONTINUE 

00  6900  1=1,3 
A I N T 1 (1) =0.0 
AINT2C1) =0.0 
OU  6600  IT=2,NZWGL 
IMUlT-1 

0T=TW6L(IT)-TWGL(IH1) 

A INTI  (IT)  = AINTl  CIMl)+0.S*OTit  (IOXCI.IT  ) + IOW(  I,  IMl)  ) 
AINT2(IT)  =AIVT2(IMl)+0.5*OT*(ZWGLS(I*tt.  IT)-t-ZWGLSC 
tl^.IMl)) 

6600  continue 

Cl=-C  (XU(I)-XP(n  ) /TOOP(I)  )*AINTt  ( N Z w G L ) 

C2  = Cl ,0/TOOP (I) ) *AtNT2(NZWGLl 
X8AR(I+1)=XICCI>1)+C1+C2 
DO  6700  IT=1,NZWGL 

Tl=-Cl-( (XOCI)-XPC!) )/TOnp(I))*AlNTl  (IT) 

T2=-C2+C1 .0/TDOP (II ) *AINT2(IT) 

XXGLS  CI  + 1 , IT)  = T 1 -f T2 
6700  CONTINUE 

Cl»  (KF(I)  / (TF(I)  *TF(I)  ) ) *AINT2(NZWGL) 

XQAR (1+4) =XIC(I>a) +C1 
DO  6600  IT= 1 , NZWGL 

Ti  = -C1*C<F(I)/(TF(I)*TF(I))  )*AINT2(IT) 

XWGLS (I+a, IT) = T 1 
6600  CONTINUE 

6900  CONTINUE 

00  7000  1=1,7 
TYPE  J»,X8AR(I) 

7000  continue 

PAUSE  »COR I C * 

00  7100  1=1,13 
ZdAR(I)=ZSA«SCI,l) 

7100  CONTINUE 

PRMT(l) =0.0 
PR*T (2) sTMAX 
PRMTC3) aHSLOW 
PRM  T (9) =0.001 
OU  7200  1*1,7 
XDOTCI) »EWS 
7200  CONTINUE 

ICRCT=0 

00  7220  I«1,NX8R 
TqOT (I) =TBAR  (I) 

7220  CONTINUE 

NM 1 =NX3R- 1 

00  7260  1*1, NMl 

OT sT8AP ( I ♦ 1 ) -T8AR  (T) 

00  7240  J*  1 , 1 3 
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zoot  u . i)  * czpars  cj.  i + n -zbars  c j,  n ) / ot 
72aj  CONTINUE 

726:j  CONTINUE 

OU  7280  1 = 1#  13 

ZoOT (I ,N*8RJ ■Z00TCI,NM1) 

7 2*0  CONTINUE 

OU  729*  1=1,13 
TYPE  * i ZOOT  ( I f 1 ) 

72<»a  CONTINUE 

PAUSE  'ZOUT' 

N U 0 T * N X 5 R 
N X B K = l 

CALL  RKGS  (PRNT,  XPaR,  XLlHT,  7,  IHLF,  FxPaR,  OXflAR,  AuX) 

N^FrsN^OP-I 

00  7400  IT*i,MZwGL 

CALL  INTPQL  ( TIhGL  f IT)  ,TBAR,NX9R,JF,TINn 
OU  7300  1*1,7 

XC(I,IT) -TINT# (XB APS  Cl. JPt l ) -X8AR3CI.jp) )*X8A«8( I, JF1 
l*XWGLSCI,IT) 

7300  CONTINUE 

7400  CONTINUE 

MXCaMZWGL 
00  7520  I a 1 # N Z W G L 
TxC(I)aTwGLCI) 

7 5 0 0 CONTINUE 

PKMT(n=0.0 
PRMT  T21 sTMAX 
PR M T (3)  aWFAST 
PRMT  (4)  *0  , 00 l 
00  7600  1*1,13 
XijOTCDaEViF 

Zwglci) »xicct*7)-zbarsci.  n 

7003  CONTINUE 

NZwGLsi 
I CRC  T = 1 

CALL  RKGS CRRMT ,ZwGL, XOOT, 1 3, IHLF.FZwGL.QZwGL. AUX) 

nZwgl=nZWGL-1 

typ£  7702 

77?a  format c * plot  corrected  results?  *,$i 

7803  ACCEPT  l , I A NS 

IFCIANS-*Y#)  7920,8100, 7900 
7903  IF(IANS-;N;)  3030,8203,8020 

8003  TYPE  1600, IANS 

GO  TO  7802 

8103  CALL  GRAPhS(NXCT,NXC,NXBR,NZNGL,  10 AUO , X S , TS , XC , TXC , 
IZ8ARS,  T8AR,ZwGLS,TiNGL) 

8203  TYp£  3300 

8303  FORMAT  ( * 00  MOPE  CORRECTIONS?  ;,S) 

8400  ACCEPT  i , IANS 

IFCIAN$-'Y')  8500,4303,8530 
8503  IF (lANS-'N')  8630,8700,8630 


8603  TYPE  1622, IAN3 
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FORMAT ( * RUN  ANOTHER  CASE?  ',$1 
3900  ACCEPT  i.IANS 

IF (lANS-'Y'l  9 000,100,9004) 

9000  IF  Cl  ANS-*N;)  9100,  <*2*0, 9100 

9100  TYPE  1300, IANS 

GO  TO  3900 
9200  STOP 

END 

l******»*#<Mr*6LQCK  DATA  PAGE  2************* 

QLOCK  OAtA 

REAL  Ka(3) .KEC3] ,KF(3) 

Dimension  XQ(3)#*Qr31f XPC3) .T00PC31 |TQ0P(3) ,TaC3) , T£(33 , 
1TF(3) ,H(3l ,0(3) 

CUMM0N/QLKl/Y0,Xr5,YP 

CUMM0N/8LK2/TD0P, TQQP 

COMMON /8LK 3/ T A , TE , TF 

CQMrtClN/gLKA/H.O 

COMnON /0LK5/ AS  AT, 6SAT 

CQMM0N/aLK6/KA,KE,Kp 

OATA  (XOCI) ,1=1, 3) >0.6, 0.3953, 0,9/ 

OATA  (XQfl) ,1= 1 ,33 /0. 58, 0.3645,0.85/ 

OAtA  (XPCH , 1*1, 3)/0.0608, 0.1 198,0.1813/ 

DATA  (TOOPfll  ,1  = 1 ,3)  /«. 0,6. 0,5’. 3/ 

OATA  (TQOP  (I)  , 1 = 1 ,3)  /0 . 25 , 0 . 54 . 0 . 65/ 

OATA  (TACX3, 1*1,33/3*0. 06/ 

DATA  CTECI) ,1=1,31/3*0.5/ 

DATA  CTFfI), 1=1, 33/3*1.0/ 

OATA  (KA(I3, 1*1, 33/3*25.0/ 

OaTa  (KE(I3, 1*1, 33>3*-0, 0045/ 

DATA  CKFCI3, t=l,3)/3*0. 16/ 

DATA  (H (I), 1*1, 3)/ 23. 64, 6. a, 3, qi/ 

OATA  (D(t) .1=1.31/9.6,2.5,1.0/ 

OATA  aSaT/0. 001123/, BSAT/0.3043/ 

ENO 

1 «****«******FUNCTiaN  S£  PAGE  3 *********** * 

Function  se(argi 

common /gLK 5/ ASAT.BSAT 

SE=ASAT,ExP(RSAT*ARG3 

RETURN 

EnO 

1****»*****«suproutine  fxct  page  a«***«***** 
subroutine  fxctctime.x.xoot) 

REAL  KA(3) ,KE(3) ,KF(31 
REAL  10(3) ,10(3) 

DIMENSION  X(20l,XOOTC203  »*0(33 .XQ(33 ,XP<33  #TOQP<33  »TDOP 
1 (3) ,Ta (3) ,TE(3) ,TF (33  »H(33  #0(33  »Y (3,33 , Th (3*33 ,VREF (33 , 
2PIN(3) ,OELR(33 ,VLT(3) 

COMMON /BLK 1/XC . XQ, XP 

C0MM0N/9LK2/T00P.T00P 

C0MM0N/8LK3/TA,TE,TF 

CUMMON/QLKfl/H.O 

C0MM0N/8LK6/K A , KE , XF 

COMMON/dLKT/Y,TH 

C0MM0N/BLK8/8H.SH 
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CQHmON/3LK9/VRFF,PIN 
COM.SOn/QLK  t l /VLT 
dei? ( i ) *3.? 

nELW(2)sBH*X(S)/H(1)*Ht3)*X(l0)/SH 
DEL* (3) B*ri*X  C8)/u(l)-H(2) *XC101/SH 
OMl*x (1 ) -$M*X  (9) /h  ( 1 ) 

OMjaX (i ) *X(91 *x  Cl  1 ) /5H 
0M3aX(l)*XC91-HC2)fiX(tt)/3H 
HO  2Cy  1*1.3 

luriiM,’1 

10  ( II  *51.3 
OU  1?>0  J * 1 . 3 

ANGsTm  cl . J) ♦HELM ( J) -HELP  (I) 

CaNGsCOS  (AMG) 

3 A '*  G a S I N ( A N G ) 

10(1) a 10 (15  *YCT, J) * (X (J*17) *CANG-X (J*i) *SANq) 

iy (i) *io(i)*y cr, j) *(X{j*t ) »cang*xcj*i7) *sangj 

10?  CONTINUE 

2u«i  Continue 

P.iHIfJaBlN  (1)  /0Ml*PlN(21  /0M2*PIN(3)  /Q*l 
OAMar>(t)*(GMt-l,0)*D(2)*(Q^2-l.  3)*O(3)*(OM3-l,0) 

Pw PuUTsO.3 
OU  300  1*1.3 

PWPQUT aP^ POUT^X  (IM  ) *10(1]  + X ( 1+ 1 7 ) * 10  ( T ) 

30?  CONTINUE 

XuOT (l) 3 (2. 5 /an} * (B^BlN-nAN-PNPOUTl 
OQ  40Q  1=1,3 

XOQT(I*x)«(r.O/TDOP(I))»(-XCl^l)-(XQCI)-*P(Il)*IO(I) »<tl 
1 + tp) 

ao?  continue 

00  530  1*1.3 

X00T(I*X|}*(l.a/TF(I))#(-X(I»4J*KFCI)*X(I*lt  J/TF(IJ) 

500  CONTINUE 

XOCT(3}*377.3*XC<n 

XOOT(9)*(0.5/6H) * (-PIN (l)/OMltO(l )*(0M1 -1.0) +XC2) *10(1)+ 
IX ( 18) *10(1 )*(H(1) /SHI* (PIN (2 J/0M2-0 (21* (0*2- 1.0) -X(3)* 
210 (2) -X (1 9) *IO(2)+PlNC3)/OM3-O(3)*(OM3-1.0)-X(a)*Tg(3) 
3-X (20) *10(3)) ) 

XQOT (10) *377. 0*X  fl  l ) 

XOOT(II) *(0.5/H(2) ) *(PIN(2) /0M2-0 (2) * (0M2-t .0) -X (3) *XQ 

1 C2J-XC19) *10(2))-.  (0. 5 /H  (3) ) *(PIN(3) /0*3-0(3) *(0*3- 1.3 
2) -X  (a)  *10  (3) -X  (20)  *10  (3)) 

00  600  1*1,3 

xooTci  + m *(i.0/TE(i))*(-(KE(i)t.sE(xa*in ) )*xci*m* 

IX (1*19) ) 

600  continue 

00  700  1*1,3 
Al*X(Ul7)*XP(I)*lQ(I) 

A2*X(I  + l)-XP(n  *10(1) 


VL7(I) *SORT (A l*AltA2*A2) 

XQOT  C I*l«)*(  l.  0/TA  CD)  * (K  A (I)  »(X  (IM)-KF(  I)  *X  (!♦!!)/ 
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XD0T(I*l7)*(i  .*/T00P  (I ) )*(-X(I*l7)«-(*Q(I)-*P(I))*ia(I)) 

au«  continue 

PETuPN 

End 

1 *«*********SU0ROUTINE  0 X C T PARE  5*********** 

SUBROUTINE  0»CT(TI^£,  X,  xOf]T,  I H|_F  , NO  I M , °«  » T ) 

DIMENSION  X(2Vi)  ,X0nT(2*)  ,PRmT(5)  ,xS(23.5*i)  ,TS(53t)  , 
l VL  T ( 3 5 

COMMON  / au  K l <)/*  S,  TS  , NXCT 
C0MM0N/8LK1 l/VLT 
00  1*0  1*1.2* 
xS(I,nxC7)*X(I) 
tJ*  CONTINUE 

Do  2*0  1*1.3 
XS(I  + 2*,  NXCT)  «VLT(I)- 
20*  CONTINUE 

TS (NXCT) *TImE 

IF(NXCT.r.E.52t)  P»MTC5)*1.* 

NXCT*NxCT*l 

seTupm 
EuD 

!•*«**** #**«SUPPnUTiNE  SLVFST  bare  6*********** 

subroutine  sLvPST(xaxR,zBAP,riM£i 

01  MENS  TON  H(3) ,0(3) . <8  AP ( 7 ) ,ZB*R(13) .T(3,3) »TH(3*3) i 
1TA(3),TEC3).TFC3).XO(3),XR(3),XP(3).OEL0(3).A(5,S). 

20(5) , WKAPEA(5) » V (33 .PIN  (3) ,VWEF(3) , ZOO T ( 1 3 , 5a 1) .TOQT( 

35*1) ,0QT(131 . TQOP ( 3 ) ,TQGP(3) 

BEAL  KA(3)  ,KE(3) ,<E(3) , 10(3) .10(3) 

C0MMQN/8LK  t/XO, XO, XP 
COMMON/8LK2/TOOP, TOqP 
COMMON/0LK3/TA.TE. TF 
C0MM0N/bLK4/H, 0 
COMMON /HLK 5 /AS  AT, 08  AT 
CUMM0N/BLK6/I<A,KE,KF 
CUMMON/dLKT/Y , TH 
CUMM0N/6LK )/»m , S M 
COMNOM/0LK9/VPEF,PIN 
CUMMON/0LK 15/ICPCT 
CQMMON/9LKt7/ZOOT,TOOT,NOOT 
C l * C 1 . WfM  c 1 3 /3H) *SIN(TH(1#2)1 
C2*(l ,0-H(l) /SH)*C0S(TH(1 ,2) ) 

THia*y  (i  .a) *S0PT(C1 *Cl*C2*C2) 

PSIl2*ATAN2(Cl ,C2) 

C1«(1.0*H(1) /SM) *SIN(TH(1,3) ) 

C2«  (l.O-M(l) /SH) *C0S(TH(1 ,3) ) 

7H13*Y(1 ,3)*S0RT(C1»C1*C2*C2) 

PSI13*ATAn2(C1,C2) 

Cl*(l.a/M(2)*l.a/H(3) )*SIN(TH(2,3)) 

C2*(-l,0/H(2)+l.*/H(3)) *COS (TH (2 . 3) ) 
YH23*0,5*Y(2.3)*SQ°T(Ci*Cl*C2*C2) 

PSl23*AT4N2(Cl.C2) 

00  22  1*1,13 
DOT  (I) *0.0 
23  CONTINUE 
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IFUCRCT)  00,8n,a,j 

4.)  CaLU  INr?QL  (TIM£,  TOQT.sjorJT,  JF.TJMT) 

00  &0  1*1,13 
OCJT  ( I j s^OQT  £ I , JF ) 

8«5  CONTINUE 

fld  z0  AR  C23 sOQTCl 3/377.3 

ZeAR(a)*0CT(3)  / 3 7 7 . 0 

OMisxeAsm-SH*zeA«(2)/H(n 

0M2*XrjAR  m *Z0AR  (21  *H  (3)  *Z0A»  (4)  /SH 
0M3*XBA»m  ♦ZdArt  (2)-H(2)*ZPAR(4)/SH 
00  14^)0  I TER* 1 ,50 
ICONV* i 

00  200  1*1.5 
0(1) *0.0 

00  l?H  J * l » 5 

A (I, 

100  CONTINUE 

2U0  CONTINUE 

ANGl*aw»Z^AR(l)/HCl)+H(3) *Z0AR(35 /SH 
ANG2s0H,ZBAfi (1) /H(l) -H(2)  * Z 3 A R (3)  / SH 
GANi*PSH2*4NGl 
GAR2*FSI13*AMG2 

Tl*(H(l)/SH)*(PlN(2) /0^2»PIN(31 /0M3J-PIN(D /0M1 
Ta«o(i)«(OMi-i‘.ai-(H(n/SH)  * cocai  *(0M2-i.8)+0(3)  *com3 
1-1.0)) 

Ti*(XBAR(2) *XBAR (2)+Z»AR(l t ) *ZBAR£l 1) ) »Y(1, l) *CQS 

1 (TmCI.  1))-(H(1)  /SH)  <t  ((x0AR(3)  *X0AR(3)  *20aR(12)  *Z0A 
2R(12))*Y(2,2)*COS(TH(2,2))t(X0AR(4)*X0AR(fl)*Z0AR(l3 
3)*ZBaO(131)*Y(3,31 *CQS (TH (3  # 3) ) ) 

Pl*X0AR(3)*X0AR(3)+Z,?AR(tn*Z9AR(12) 

P2*XPAR  (?)  *20AB(i2i  -XBAR  (3)  *20 AR  CU) 

P 3*X0APC2) *X0AR(4)*Z0AR(ti)*Z3AR£l3) 

P43XBAR (?) *Z0AR(i3)-X8AR (a) *29 AR (1 U 
P5*X0aR(3) *Z9AO(13)-X8AR(4) #ZaAR(l2) 
Pb*X8AR(3)*X0AR(4)*Z0AR(l2)«ZaAR((3) 

T4*THi2*(Pl*C0S (GAMi) tP2*SlN(GAMl) ) 
T5*YNl3*(P3*C0S(GA'‘21*p4*SlN(GAN2)) 

Tb*2.0*(H(l)/SH) *y (2,3)*C0S(Th(2,3) )*(P5*SIN(Z3AR(3)) 
l-Po*COS(Z0AR(3))) 

8(1)  *-(Tl*Ta*T3*T4*T5*-Tb)  +2.0*0H*OOT(2) 
Tl«2,0*ZBAR(U)*y(t,  i)*C0S(TH(t,  l)) 

T2*yPl2*f20AP (12) »C0S(GAMi)-X8AR (3) *SIN(GaM1)) 
T3«YH13*(Z0AR(13) *C0S(GAM2)-XBAR(4) *SIN(GA*2) ) 
A(l,l)«TUT2*T3 

Tl«-2,0*(H(l)/SH)*Z8AR(l2)*Y(2,2) *CGS  (TH  (2 , 2 ) ) 
T2*yH(2.(Z0AR (l 1) *COS (GAMi) +X8AR(2) «SIN(GAMl) ) 
T3*-2.0*(H(1)/SH) *Y(2,3)»C0S(TH(2,3))»(X8aR(4)*SIN( 
1Z8AR(3))*Z8AR(13)*cOS(20AR(3))) 

A ( 1 # 2) *T1+T2  + T3 

Ti*-2.0*(H(1)/SH) aZ8A«(13) *Y(3»3)*C0S(TH(3,3)) 
T2*YH13*  (Z0AR(1  n*CaS(GA*2)  4-X3AR(2)  *SIN(GAH2)  ) 
T3*2.3*(M(1)/SH) *Y(2,3) *C0S(TH(2,3) )*(X3AR(3)«SIN( 
128a«(3) )-Z8AR(l2) *COS(2BaP£3) ) ) 

* C l » 3) *T l ♦T2  + T3 

li 
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Tl*  (YH12*BH/H  (1 ) ) * r-pi*SlN(GAMl)  ^oCOS  (GA*U) ) 

T2s  CVH13*i«/H(in  *C-P3*SI*(GAM2)  tP«*C0S  (GA*2) ) 
*(l»fl)aTltT3 

Tl  = (YHl2*Hf35/SHj  *C-Pt*SlN(GAMl)*P2<rCaS(GA*U  ) 

T2* (YHl3*Hf2i /SH) * fP3«SIN(GAM2)-Pa*C0SfG4M3) ) 
T3=(2.^*Hcn/SH3*Yf2,3UC0S(THf2,3n*(p5*C03(ZBAP(3 

1 )  ) +P&*STNfZ8AP  (3) ) ) 

A(t,5)*n*T2*T3 
GA^isTH(2, l)-ANGt 
Gam2sTh(3, 1 ) -ANG2 
GAMjsP5I23*Z8AP  (3) 

Tl*i).5*(pIN(P)  / (0M2*H(2))-PINC1]  /(0^3*H(3)  3 ) 
T2a2.S*(-Of2)*(aM2-l.ai/H(2);of3)*fOP3-l.3l /H(3) 3 
T3*0.5* C-(**ar  (3)  *Y8Af?(3)*ZSAP(t2)*ZBAP(12)  )*Y(2,2)  • 
lC05(TH(2,?n/MC2)  ♦(Y8AS  ( a ) * X a AR  ( a ) ♦ Z3  aP  ( 1 3 ) *ZB  A*  C t 3 3 

2) *Y(3,3!*CUS(THf3.3))/H(3)) 

T4»-(3.5*V  (2. 11  / H(8) ) * (Pt  *C0S  (GAMn-P2*STN(GAMin 

T5«C9.5*V(3,n/Hf3n*  (P3*C0SfGAM2)-pa*siN(GAY'2n 
T<>«YH23,  (Pb»C03  (GAH3l«P5*siN(GAM3) ) 

B(2)a-(TUT2>T3  + Ta>T5>Tfe)  tOQT(4) 

Tl*-a,5i»(YC2.  n/H(2))*(ZBAR(l21*C0S(GAMn+X8AR(31  *SIN( 
l G A M l ) ) 

T2av).S*(Y(3, 1)  /H(3n  *(Z9ARC13)  aCOS  CGAM2]  ♦XBaR  (4  ] *S  I M ( 
1GAM2] } 
a(2,i3sTi-»-t2 

Tls-ZuARCiP) * Y (2*21 *C0S(TH(2.21 ) /H(2) 

T2«-(0.5*Y(2,  n/H(2))*(ZBAR(intC0S(GAMn-xaAR(2)  * S I N ( 

1 G A M l ) ) 

T3«YM23*(ZPAP(t3}*C0S(G4H3)*XHAR(<O*SIN(GAM3)) 

A (2 , 2) s T i +T2*T3 

Tl«ZBAM(l3)*Y(3,3) *C0S(TH(3,3) 1 /H(3) 

T2=(0.5*vC3,n/H(31  ) *(ZBAR(U  J*COSCGAH2)-XBaR(2J  *SrN( 
IGAM2) ) 

T3aYH23*fZBAP(l2) *C0S(GAM3)-X8AR(31  *3IM(GAH3) ) 

A (2,33  *T1+T2>T3 

Tl*(0,5*Y(2,l)*eH/(H(23*H(l)))#(-Pt*SIVCGAMi)-P2*COS( 
IG**1!)] 

T2*-C0.S*Y(3. 1) *BH/(H(3) *HC13) ) * (-P3 *3 t N (G AM2) -pa*COS ( 
IGAM2))  | 

A(2,ft) aTl*T2 

Tl*  (0.5*V(2, 1) *H(3)/  (H(2) *SH] )*(-Pl *3I^(GA^i)-P2*C0S( 
1GAM1) ) 

T2*(0.5*Y(3, 1) *H(2) /(H(3)*SH])*(-P3*SI\(GAM2)-P4*C0S( 
1GAH2] ] 

T3*rH23*(-Pb*SIN(GAM33-P5*C0S(GAM33 ) 

A(2,5)*TUT2  + T3 
0£LR  (1) 12,0 
DEL*  C23  »ANGl 
OElR  (33  sANG2 
00  4910  1*1,3 
:q(U»0.« 

00  300  J>1.3 

ANGsTHfl, J)>nELR(J)-OELR(I) 

CA^GaCOSfA^G) 


• f 


J 


J 
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320 
420 

*ae 


620 

720 

f|J0 

900 


J 


9b3 


1102 

1202 


1302 


1403 

1502 

1602 


i 


SaNGsSInCAMO) 

I<3£l)*ry(I)  ♦’Yfl,  J)  *£XBAR(Jt-n  *CANG«-Z049£J*10)  *SaNG) 

CU^TINUc 
CONTINUE 
OU  502  I ■ 1 . 3 

8 (I*2)  mZBA»  Ct^iai  - f XQ(I) -XP(IJ 1 *IQ (II ♦T30P  (I) *OOT  (1  + 101 
CONTINUE 
Oo  902  1*1.3 
00  202  J * 1 . 3 

ANG*Trt(I. J)*0£LP£J)-OElRCn 
IF(I-J)  70fl, *,0", 700 

Act+2,  Jj«-f  l.a-cxorn-xp(ii)  *yri.n  *sim(thcii  in  i 

00  TO  220 

A (1+2,  J)  »(XQf  Il-XP(l))  *Y(I,  J)  •SIN(ANO) 

CONTINUE 
CoNT INUE 

GA*l*TM(l  ,2)  +ANG1 
Ga“2*Th(1,3) ♦AnG? 

Tl*r  (1 ,2)*(-xaAR(3)  ■SIN(GAMH+zaAR(l2)  aCOS(GAMI)) 
T2iY£l,31*(-X8A* (41 »SIN(GAH21 +Z3AH(13) *C0S(GAH2) ) 

A (3, 4]  ■(XQm-XPfll  t(Ti*Tai 

Ti*H(3)»T(l,2)*C-X«A0(3)*SlN(GAHl)+Z9A9£l2)*COSCGAMin 

T2*H (2) *Y(l#3)*C-X«AP(4)*SlN(GAM2)^Z9A9(t31*C0S£GAM2)) 

A£3.5J*£(XQ(1)-XP£1) 1/SHl *(T1-T21 

GA^l*TH£l,21-ANGl 

GAm2*Th£2.3)-29AP£3] 

A(a#4J*-(XQ(3)-XPC2))*Yf3,l)+£9H/H(l))*£-X0AP(2)*SlN£GAMt)+ 
lZ9AP(ll) *C0S£GAM11 J 

Tl«Y(2.n*(H(31/SHl«(..X9AR(21*SlN£GAMl)i-Z9A9(in*caS(GA*t)l 
T2«Y(2.31  *f-X8A»C4)»SlN(GAM2)+Z8AR(l3) •COSCGAna) ) 
A£4,5)«-(X0(2)-XP£2))*(TWT2) 

GANi«TM(l,3)-ANG? 

GA02*TH£2, 3) *Z0AR  (3) 

A£5,4)«-(X0(3)-XP£3))*Y(1#3]*£9H/H(ll)*(-X8AR£3)*SIN(GAM1) 

l+Z0AR£lll+COS(GAMt)) 

Tl«£Y(ll3)*H£2)/SNl*(-Y9AR£21*SIN£GAMn4-Z0AR(in  *COS  (G  A H 1 ) 1 
T2*Y(2,31  *(-xe*»(31  t3lN(GAH2l  +Z9AR(13J  *C0S(GAM2n 
A£5|5)t£XQ(31*XP£31}*£Tl»T2) 

CALL  LE3TiF£A, 1 ,5,5,a.0,WKAREA.  IER1 
DO  1230  1*1,5 

IF £A8S(8£ I)) -2.0001)  1200,1200,1130 

icoNVig 

CONTINUE 

Oo  1320  1*1,3 

Z8AR(I*l0) *8£I)+Z8AR(I*10) 

CONTINUE 

ZB AR ( i ) *Z8AR(i)+B  £4] 

Z3AR(3) *Z8 ar (3) +g (5) 

IF(ICONV)  1 402 . 1400 , 1600 

CONTINUE 

TYPE  1500 

FORMAT £*  N0N-C0NVEPGENCE-ANGLE3  ANO  EOP'l 
STOP 

0EL2 (1) *0.0 


md 


r . — — 3 
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OELR (2) *art*Z94  0 ( l ) /H(n*H(3)  *ZRAR(3)/SH 
0tt-R(3)sau*zaA»  Cl ) 'H(t )-H(2) *Z9aR(3)/SW 
00  180*  1st, 3 
I0CI3  *oi.^ 

IQCI) *3.0 
00  1 7 J»1  ,3 

Ai,GaTM(I,J)*0£LP(J)-OELR(I) 

CAVGsCOSf AMG1 
SA*‘G*SIN  C AnG) 

io(n*ro(n+Y(i#j)*(Z0ARCj+i*)  *cang-xsar  c j ♦ i ) *sang) 

IQ(I)*!Q(I)+Y(I,  *CANG*Z9a»CJ*U)  -SANG) 

17«*|  CONTINUE 

l<33j  CgMflNijE 

00  IR03  I st,  3 
VO*Z0A»CI*l 3) *XP(I) -TQfl) 
vu«xeAO(i»ij-xP(ij*io(i) 

V (I  ) *SORT  (Vi)*VO  + VQ*VO) 

IR3J  CONTINUE 

00  2503  I s 1 , 3 

00  2203  ITER*l,53 
ICONVst 

ANUM*(<E(n  t(KACI)  *KF(I)/TF(m+SE(Z0AS(I  + 4))  ) *Z9aR 

1 (I*4)-.<A  (I)  *(XaA*(I*4)-V(I)  *VR£F(I)  ) *T£(I)  *00T(I*4)  *TA( 
ID  -Oqt  (1*7) 

OENOMa-C<Ef I) *(KA (I) *KF(I)/TF(I))+f3SAT*Z9AR(I+4) *1,0 
l)*SE(Z3AR(I*4))) 

OlTaANUN/OENOM 

1 IF  (A 05  (OID  -3.3051 1)  2100,2103,2000 

2 I C 0 N V * 3 

2130  ZBAR (I*  i) sZ3AR  (1*4) *0IT 
IF (ICON  7)2? 3 3, 2230, 240 3 
2230  CONTINUE 

TYPE  2 3J0,  I 

233.)  FORMATC'  CONVERGENCE  FAILUPE-EFO  ;,I1) 

STOP 

2430  Z3AR  CI+7) =X4 (I) - (XBAP  (I*4)-KF(I) -Z9AR ( I +4 ) /TF ( I ) -V  ( I ) 
1*VREF(I) )-TA(I) -OOTCI+7) 

2530  CONTINUE 

2630  RETURN 

ENO 

J ******** ***SURROUT  IN£  FX 9 A R PAGE  7********** 

SUBROUTINE  FX8AR (TIME, X0AR, xDQT) 

PE At  KA(3) ,KE(3) ,KF(3) ,10(3) ,IQ(3) 

OINENSION  X0AR(7) ,X00T(7) , Z9 AR ( 1 3) , XO (3) , XO (3) , XP (3) , 
ITD0PC3) ,T00P(3) ,TA(3) , TE  ( 3 ) . TF  ( 3 ) . H ( 3)  , 0 (3 ) , Y C 3 » 3 ) , 
2TH(3,3) ,VP£F(3) ,PTN(3) ,0ELR(3) 

CQMMON/8LK1/XO, XQ, XP 
COMMON/9LK2/TOOP,TOOP 
COMMON/ 9LK3/T A, TE,TF 
COMMON /8Li< 4 /H,  (5 
COMMON/0LK6/KA,KE,KF 
CQMH0N/oLK7/Y,TH 
CQMmQn/3LK8/«h, 3W 
COMMON /0LK 9 /V REF, PIN 
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C{J**0n/3L'<127Z9AR 

CALL  SIVFST(X9A*,Z9aR,TIME} 

0£LH ( 1 ) *3 

OELR (2)s8H#ZBARCt) /HCi)+W(3) *Z^AR(3) /SH 
DEL*  (31*8H*ZaAS(l)/Htn-Mc2)*ZaARC3)/5H 
0MiaxeARfn-3w*Z8A3(2]/Htl) 

OH2«X8A«n)»Z3ARf23+H(3)  *ZBAO(a)/SH 
0M3*K8A2f  l)4-ZdA«f2l-H(2)  *z9aR(4)/SH 

OU  200  I *1.3 
10(1) *0.0 
IQ  C 1 3 *0.0 
ou  j*i,3 

ANGsTHd  , J)  ♦OELR(  J)  -OELR  (I) 

CA^GaCOS  ( AKIG) 

SAMQaSIN ( AMG1 

I0U)«I0(I)+Y(I.J) • (Z8AR  C J ♦ l 3 ) *CANG-T3AR  C J * 1 3 *SANG) 

iorij«xa(n*vci.  j)*(*9AR(j+n«CANG*z9A»(j*ta)  asangi 

\ l?  CONTINUE 

2*0  C(jnTInu£ 

PWRIN«PIN(1 ) /ON  1 ♦PIN (2) /0M2tPiN(3) >0*3 
nAN.OCn*(OHl-l.ffl)*Q(2)*(QM2-l’.3)AO(3)*(OM3-t.0) 
PwROUTa'j.3 
OU  400  1*1.3 

o «tRGUT*PwR0ijT  + X8AR(I  + l ) *lQ(n*Z8AR(I*10)*IO(I) 

400  CONTINUE 

X 0 0 T (1)  *(0.5/8H)  * (PwRIN-OaM-pwPOUT) 

00  500  1*1.3 

XOOT(I  + l)sa'.0/TOQPCln*C-X3A9(IM)-(XOCn-XPCl))* 
lI0(I)+Z8A»CIt4)) 

500  CONTINUE 

00  600  1*1.3, 

X00T(I*4)  *(l.0/TF(I) ) *(-X3AR(I*a)*KF(I) *Z8AR(I*4) 

l/TFCl) ) 

600  Continue 

return 
End 

!**********$ U9RQUTINE  0x8aR  PaGE  8*******«*** 
subroutine  oxqar (time, xbar, xoot, ihlf. noim,prmt) 

01  men si on  X8AR(7) ,X00T(7) ,PRMT(S) .XBARSCZ,50t) ,Z8ARS( 
113,501) , ZB AR( 13) ,TBAR  (531) 

CQMNON/dLKte/ZBAR 

C0MM0N/8LX13/XBARS, ZBARS. T3AR, NXBR 

Do  100  1*1.7 

X8ARS (I,NXBR) aXBAR(I) 

100  continue 

00  200  1*1  , t3 

ZBARS (I,NX8R) *Z3AW  (I) 

200  CONTINUE 

TBAR (nxBR) *TIME 

IF(NXBR.GE,50n  PHMJ  (5) * l .0 

Nx8R«Sx3Rtl 

RETURN 

EnO 


l 


SUBROUTINE  FZ'AGL  page  r 
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SUShOUTINE  FZWGL(TImE,ZwGL,xPOT) 

REAL  KA(3).KF(3).*P(3),IOC3),IO(3) 

01  PENSION  Z wGL  ( 1 3 ) . X00T(13) ,X0(3) ,X0(3) ,XP(3)  , T0OP(3)  , 
ITQQP  C3).TAC3).TE(3)»TF(3),Hf3Jf0C3),Y(3.3).TH(3»3l, 
2VR£F(3).PIN(3) , OELR (31 , XB ARS (7 , 53 1) , Z* APS  ( 1 3 , 50 1) » 
3TBAp  (501 ) . XC(7,5«J1)  , TXC  (501)  ,X  (7)  , Z(13) 

C0MM0N/8LK l /XO  , XQ , XP 

cuMMON/ai.'xa/Tuop.  toop 
CQM10M/BLK3/TA, TE, TP 
CUMM0N/8Li<a/H#0 
CCJMM0N/9LK6/KA.KF.KF 
CUMH0N/3LK7/Y, TH 
CaMH0N/9LK3/«H, SM 
COMMON /8LK0/VREF, PIN 
CO-lHON/aLX  13/X8APS,  Z3ARS.TBAR,NXflR 
COMMON/aLKU/XC,  TXC,NXC 
C0MM0N/8LK15/ICRCT 
IFCICRCH  130,100,300 

CALL  I'iTPOL  (T1M£,  TBAP.NXaR.JFINO,  TINT) 

00  200  1*1,7 

X(I)aTlNT*fXBA9Sri,jFlNO<t>l)-X8ARSCI,JFTNOn*XaARS( 

1 I*  JFINO) 

continue 

GU  TO  500 

call  intpolctime, txc.nxc, jfino .TINT) 

00  400  1=1,7 

X(Ij*riNT*(XC(I,  JFtMO+n«XCCt,  JFINO)  1+XCCI,  JFXN01 
CONTINUE 

CALL  lNTpOL(TlME,TBAR,NxaR, JFImQ, TINT) 

00  600  1*1,13 

Z(t) ■TINT*(Z0A»3(I, jFlN0*l)-Z8APS(t, JFINO) )»ZBARS( 

1 If JFIN0)+ZWSL(I) 

continue 

OELN  (1) *0,0 

06LW(2)*9h#Z(1)/H(1)*m(3)»Z(3)/SH 
OELR (3) sBh*Z  fl) /H(l)-H(2J *Z(3)  /S* 

0m1sXC1)-SH*Z(2)/H(1) 

0M2*x  (l)>Z(25+H(3)*z(«)/SH 
OP3sX(l)t.Z(23-H(2)*Z(4)/SH 
00  903  1*1,3 
10 (I) *0,0 
IQ(I)*3.0 
00  730  J * 1 , 3 

ANS«TH(I, J) *OEL«(J) -OELR  (I) 

CANG=COS (ANG) 

SANGsSlN(ANG) 

10(1) *10(1) ♦Y(I, J) * (Z (J*l0) oCANG-X (J*l) *Sang) 
XQ(I)*IQ(I)+Y(I,J)*(X(J*l]*CANG+Z(J+t3)*SANG) 

CONTINUE 

Continue 

XOQT(n  «377.0*ZWGL(2) 

XuOT(2)  s(S'.5/3H)*C-PlN(l)/0Mt+0(t)*(0, -11-1.0 )4-X  (2)* 
iIO(i)+Z(tl)*IO(l)*(H(l)/3H)*(PIN(2)/QH2-0(2) *(0M2 
2-1 .0)-X(3) *IQ(2)-Z(12) *10 (2) ♦PIN (3) /qm3-0(3) *(0*3 
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3-t,P)-X(4)  *I«J(3)-Z(13)*IQ(31)) 

X00T (31 a377,0*ZnGL fa) 

XQPTCa)  = (3.5/h(2)  ) * C p I n*  (2)  /0M2-DC2)  »(0M2-1 .0)-x  (3)  * 
in(2)-Z(l2)*X0(2))-(3.5/H(3))*(PlN(3)  /0M3-0  (3)  * C 
20H3-1 .2) -x (a) *10 (3J -z ( 13) *10(3)  ) 

00  900  1=1.3 

<oor(i>4)  = a.53/TE(i))*(-(Kecn>3E(Z(r  + a)n*zci+^i* 
lZ(I+7)) 

CONTINUE 

no  ic?aa  1*1,3 

*i«z(i  + ta)*xP(n*io(i) 

42*X(I+1)-XP(I) *10(1) 

VtTsSUG  r C A 1 *Al+A2*A2) 

X00T(I  + 7)*f  l.a/T»  (I)  ) *(K4(()  *(X  (Itaj-XFd)  *Z(I  + 4) 

l/TF(X)-VLT*VRCF(I))-Z(l+7)) 

CON  T iNlJg 
OQ  113?  1=1,1 

X&OT(I*10)»(l .«/TQOP(I) )*(-Z(I*10)*(XQ(I)-XP(I))*IQ 

1(D) 

CQN  r INUE 
e>E  TljPN 
ETiO 

J ••••••••••SUBROUTINE  OZwCL  PAGE  13********** 

SUR ROUT  I Mg  OZWGU  (TIME , ZWGL, XOOT, IHLF.NOIM.PRMT) 

Dimension  ZWGL(13)  ,XD0T(13)  ,PPMT(5)  ,ZwGLS(l3»S0n  , 
lTwGL(Sai) 

COMMON/SLK »6/ZWGLS,TNGL,NZWGL 
DO  100  1=1,13 

ZrfGLS(I,NZW<SL)«ZwGL(D 

CONTINUE 

TrtGt  (NZWGU) sTIME 

IF (NZwGL.G6.501)  PPMT(5)=1,0 

NZNGL*NZWGL+t 

RETURN 

End 

•**********suqpqutine  intpoo  page  ti*********** 

SUBROUTINE  INTPOL  (TIME.T.N, JF I NO, TINT) 

DIMENSION  T(N) 

OU  100  JTs 1 , N 
JFINOsJT 

IF (A8S  CTIME-T(JT) ) -1 .3E-05)  230,200,10* 

CONTINUE 
GO  TO  33? 

TINT* a,0 
RETURN 
iNM  1 aN-l 

00  600  JTal.NMl 
JFINOsJT 

IFCTIME-T(JT))  600,633,400 
IF  (TIME-T(JTM)  }500,600, 600 

TNUM.TIME-T CJT) 

TOENOM=T(JT+l)-T(JT) 

TINTsTNUM/TOENOH 
P £ T u P N 
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1300 


1103 


CONTINUE 

OlFF  = r ImE-T  ( JFINn) 

TYPE  730. TIME, TCJFInO) ,0IFF 
FUPMAtC'  ;,3C2X,El5.8) ) 

TYPE  830 

FysNArc*  »**interpolatton  failure***') 

PAUSE  'INTPQL' 

STOP 

EnO 

J**********SU(3RQUTTnE  GRAPHS  PAGE  12********** 
SUBROUTINE  GRAPHS (NxCT.NXhAT, NZB A R, nZwGL,  IRAUD, XS, TS, 
1XHAT, TXMAT, Z8AR, TZ3AR, Z^GL, TZWGL) 

OIhensTQN  XS(23,NxCT),TSfNXCT).xHATC7,vJXHAT),rXHATC 
lNXHAT) ,Z0AO  C13,MZBaR) , TZBAP  tNZBAR) , ZWGL  C13.NZWGL) , 
2TZWGL(NZHGU  ,PLT(53l)  ,ZHAT(S0t) 

TYPE  200 

FORMAT!'  which  STATE  hIlL  RE  PLQTTEO?  ',*) 

ACCEPT  303, IPLT 
FORMATCI) 

IFCIPLT-7)  300,803,400 
IFCIPLT-23)  1300,1300,500 
00  600  I*1.NXCT 
PLTCI) =XS  CIPLT,  I) 

continue 

CAUL  INITT(IBAIJO) 
call  qinitt 
CALL  NPTS(NXCT) 

CALL  XFRMC4) 

CALL  yF*m(4) 

CALL  CHECK  (TS,PLT) 

CALL  OSPLAY(TS.PLT) 

CALL  Frame 
CALL  aNMOOE 
ACCEPT  700 , IANS 
FORMAT CAt ) 

GO  TO  1300 
AM  IN  a I , yE  30 
AMAXs-I ,0630 
00  900  I a l , NXCT 

IF(AHIN.GT.XS(IPLT, n ) AMINaXS  CIPLT,  I) 

IF  (XS  (IPLT,  I)  .GT’.AHAX)  AHAXaXS  CIPLT,  I) 

CONTINUE 

Oo  1000  Ial.NXHAT 

IFCAMIN.GT.XHATCIPLT.I))  AMIN.XHATCIPLT.I) 

IF  CxHATCIPLT,  I)  .GT'.AHAX)  AMAXsXHAT  CIPLT,  I) 

CONTINUE 

00  1100  Ial.NXCT 
DLT (I) aXS  CIPLT  , I) 

CONTINUE 

CALL  INITTCI3AU0) 

CALL  QINITT 

CALL  OLImYCAMIN,AMAX) 

CALL  NPTS(NXCT) 

CALL  XFRM(tt) 


Call  YFw'MU) 

ChI.L  CH£C:<CTS,blT3 
Call  uSFLaV (To.^LTI 
CALL  F9AM£ 

OU  1 2 >J 3 1st  .MXHAT 

°C t C 1 3 sx^ATtiPi.r,  n 
ieou  cL'*irif.UE 

call  l!nEc3*o 

Call  0LImY  (A*IN,  a,-a*) 

call  nptscmxwat) 
call  check ( txhat.plt) 

CALL  CpL-^  (T*HAT,PlT) 

call  an-iooe 
accept  7'n.[A\S 
gc  Tn 

13?,J  JPL'sIPuT-7 

00  l<At:0  1st,  MZwGl. 

Tl^EsTZ  MGLf  I) 

CALL  INT°CjL(TIME<TZiiA^,NZ3AC?<JFIM0,TlNT) 

ZMATtl)aTlMT*(ZbAKfjPLT,JFlNn+n-2RAi?(JPLT,jFIN0J) 

i •►zcak  r jplt,  jfino)  +zwgl  c Jplt.  n 

lAP/j  cu^T  rr.ufc 

a M I N = 1 . C E 3 P 
A:iAx=-i  .<*£ 3? 
no  1 5vJP  T = l , M X c T 

If  (AHXM.GT.XSCIPLT.  in  AMiNsys  CIPLT,  T) 

Ip  SCI  PL  r,n. ST ‘.AM  AX)  AHAXsXSCIPLT,  n 
istn  continue 

Ou  160P  Ist.HZHGL 
IF  (AHIN.GT.ZHAT  c rn  AMIMsZHATfn 
If(ZHATCn.GT.4MAXl  AHAX»ZHAT(I) 

it,?u  continue 

OU  170?  I=1,NXCI 
PlT  CI) s<5 CIPLT  . n 
l 7 r-* «:  CQNflNuE 

CALL  INITT(I«AMQ) 
call  aiNirr 

CALL  ULIMY (A“IN, ANA*) 

CALL  NPTS(NXCn 
CALL  XF*MC«J 
CALL  Y FK h ( a ) 

call  check its , plt) 
call  osplay cts.hlti 
call  Frame 

CALL  LINE(3«) 

call  0 L I m Y CAMIN, amax! 

Call  nptscnzwol) 
call  chec<ctz*gl.zhaT) 

Call  CPLOTCTZrtSL.ZHAt) 

call  ANf-ooE 

ACCEPT  7f*a,IANS 
1 TYPE  19<n 

19PO  FQPM  A T ( * MORE  plots?  *',sl 

24."»J  ACCEPT  7PP.IANS 
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r. 


IFdAMS-'Y'l  2100, 12*, 2103 
210a  IF (lANS-'N'J  2200, 2«00, 2230 

22«2  TYPE  2330,  T^'S 

2304  FQPNATC'  ??  * , A l , * ??  type  again  *,S) 

GO  TO  2030 
240y  R£TyPiN 

E NO 


213 


REFERENCES 

1.  E.  W.  Kimbark,  Power  System  Stability,  Vol.  1,  John  Wiley,  New  York,  1948. 

2.  J.  E.  Van  Ness,  Hewlon  Zinmer,  and  Mehmet  Cultu,  "Reduction  of  Dynamic 
Models  of  Power  Systems,"  1973  PICA  Proc. , pp.  105-112,  1973. 

3.  J.  M.  Undrill  and  A.  E.  Turner,  "Construction  of  Power  System  Electro- 
mechanical Equivalents  by  Modal  Analysis,"  IEEE  Trans.  Power  Apparatus 
and  Systems,  Vol.  PAS-90,  pp.  2049-2059,  1971, 

4.  J.  M.  Undrill,  J.  A.  Casazza,  E.  M.  Gulachenski,  and  L.  K.  Kirchmayer, 
"Electromechanical  Equivalents  for  Use  in  Power  System  Stability 
Studies,"  IEEE  Trans.  Power  Apparatus  and  Systems,  Vol.  PAS-90, 

pp.  2060-2071,  1971. 

5.  R.  W.  De  Mello,  R.  Podmore,  and  K.  N.  Stanton,  "Coherency  Based  Dynamic 
Equivalents  for  Transient  Stability  Studies,"  EPRI  Final  Report, 

1974. 

6.  J.  H.  Chow,  J.  J.  Allemong,  and  P.  V.  Kokotovic,  "Singular  Perturbation 
Analysis  of  Systems  with  Sustained  High  Frequency  Oscillations," 

Automat ica,  Vol.  14,  pp.  271-279,  1978. 

7.  R.  P.  Schulz,  A.  E.  Turner,  and  D.  N.  Ewart,  "Long  Term  Power  System 
Dynamics,"  EPRI  Final  Report,  1974. 

8.  0.  W.  Hanson,  C.  J.  Goodwin,  and  P.  L.  Dandeno,  "Influence  of  Excitation 
and  Speed  Control  Parameters  in  Stabilizing  Intersystem  Oscillations," 

IEEE  Trans.  Power  Apparatus  and  Systems,  Vol.  PAS-87,  pp.  1306-1313, 

1968. 

9.  J.  M.  Undrill,  "Equipment  and  Load  Modeling  in  Power  System  Dynamic 
Simulation,"  Systems  Engineering  for  Power:  Status  and  Prospects, 

ERDA,  pp.  394-420,  1975. 

10.  P.  V.  Kokotovic,  R.  E.  O'Malley,  Jr.,  and  P.  Sannuti,  "Singular 
Perturbations  and  Order  Reduction  in  Control  Theory  - An  Overview," 
Automatica,  Vol.  12,  pp.  123-132,  1976. 

11.  K.  W.  Chang,  "Singular  Perturbations  of  a General  Boundary  Value  Problem," 
SIAM  J.  Math.  Anal. , Vo..  3,  pp.  520-526,  1972. 

12.  P.  V.  Kokotovic,  "A  Riccati  Equation  for  Block-Diagonalization  of 
Ill-Conditioned  Systems,"  IEEE  Trans.  Automatic  Control,  Vol.  AC-20, 
pp.  312-814,  1975. 

13.  IEEE  Conmittee  Report,  "Computer  Representation  of  Excitation  Systems," 
IEEE  Trans.  Power  Apparatus  and  Systems,  Vol.  PAS-87,  pp.  1460-1464, 

1968. 


I 


214 


14.  A.  B.  Vasil'eva,  "Asymptotic  Behavior  of  Solutions  of  Certain  Problems 
for  Ordinary  Non-Linear  Differential  Equations  with  a Small  Parameter 
Multiplying  the  Highest  Derivatives."  Usp.  Mat.  Nauk. . Vol.  18, 

pp.  15-86,  1963.  

15.  K.  N.  Stanton,  "Dynamic  Energy  Balance  Studies  for  Simulation  of 
Power-Frequency  Transients,"  1971  PICA  Proc.,pp.  173-179,  1971. 

16.  E.  W.  Kimbark,  Power  System  Stability,  Vol.  3,  John  Wiley,  New  York,  1956. 

17.  IEEE  Committee  Report,  "Recommended  Phasor  Diagram  for  Synchronous 
Machines,"  IEEE  Trans.  Power  Apparatus  and  Systems,  Vol.  PAS-88, 
pp.  1593-1610,  1969. 


; 


18.  G.  Shackshaft,  "General-Purpose  Turbo-Alternator  Model,"  Proc . IEE  . 
Vol.  110,  pp.  703-713,  1963. 

19.  P.  L.  Dandeno,  P.  Kundur,  and  R.  P.  Schulz,  "Recent  Trends  and 
Progress  in  Synchronous  Machine  Modeling  in  the  Electric  Utility 
Industry,"  Proc.  IEEE,  Vol.  62,  pp.  941-950,  1974. 

20.  D.  N.  Ewart  and  R.  P.  Schulz,  "FACE  Multi-Machine  Power  System 
Simulator  Program,"  1969  PICA  Proc.,  pp.  133-153,  1969. 

21.  A.  W.  Rankin,  "Per-Unit  Impedances  of  Synchronous  Machines," 

A IEE  Trans.,  Vol.  64,  pp.  569-573,  1945. 


i 


A 


* 


215 

VITA 

John  Jay  Allemong  was  bom  in  Harvey,  Illinois  on  January  23, 
1951.  He  received  the  BSEE  degree  with  Highest  Honors  and  University 
Honors  and  the  MSEE  degree  both  from  the  University  of  Illinois,  Urbana- 
Champaign  in  June  1973  and  May  1974  respectively. 

He  was  voted  the  Eta  Kappa  Nu  Outstanding  Senior  in  Electrical 
Engineering  in  1973.  From  February  1973-February  1974  he  held  an  RCA 
Fellowship  in  Electrical  Engineering.  From  July  1974-August  1975  he  was 
employed  as  an  Electrical  Engineer  by  the  consulting  firm  of  Sargent  and 
Lundy  Engineers,  Chicago,  Illinois. 

Since  August  1975  he  has  held  the  positions  of  Instructor  in 
the  Department  of  Electrical  Engineering  and  Graduate  Research  Assistant 
in  the  Coordinated  Science  Laboratory  of  the  University  of  Illinois.  In 
April  1977  he  received  the  Harold  L.  Olesen  Award  for  excellence  in  under- 
graduate teaching  by  a graduate  student.  His  interests  are  power  systems 
and  the  application  of  computers  to  power  systems. 

Mr.  Allemong  is  a student  member  of  the  Institute  of  Electrical 
and  Electronics  Engineers. 


1 


